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ABSTRACT 



A study is inade of the thermally induced circulation in Lake Michigan. 
The surface temperature of the lake is approximated as a concave cosine func- 
tion during the warming-up seasons and as a convex cosine function during 
the cooling-off seasons. The bottom "boundary has very stable temperature all 
year round and a small heat flux is allowed to flow through the side bound- 
aries. The free surface is assumed to be free from wind stresses. 

It is noticed from a literature survey that for small: thermal Rossby num- 
ber and small Ekman number, such as found in the Lake Michigan model, flow in 
the lower symmetrical regime is expected to occur* Therefore, in this theo- 
retical study, the first step is to linearize the system of nonlinear govern- 
ing equations of the flow field by expanding all nondimensional dependent 
variables in powers of the small thermal Rossby number. Then, for every order 
of the thermal Rossby number in the expansion, a matching additive type bound- 
ary layer analysis is used to obtain inner and outer approximate solutions of 
the flow field. Under the scheme of double perturbation expansion, ignoring 
four small corners, the temperature distribution and the velocity field can be 
systematically solved for all regions. 

Two types of general thermal currents have been found from the analytic 
solutions in our theoretical study. The type-A circulation has a cyclonic 
meridional velocity and two cells of zonal and vertical flow in the cross - 
section with sinking motion at the middle and upwellings on both edges of the 
lake. This type of circulation is expected to occur during the winter-cooling 
period, usually from January to March, and during the summer -heating period, 
which is usually the period between the disappearance of the thermal bar and 
the period of strong stratification. The type-B circulation takes place gen- 
erally in late March (denoted as the spring-heating period) or late November 
and December (denoted as the autumn-cooling period). Its flow is exactly the 
reverse of the type-A circulation. 

The type-C circulation occurs during the spring and the autumn thermal 
bar periods when two weakly stratified waters, one of which is in direct strati- 
fication and above ^°C as the other is iniinverse stratification and below ^^^C, 
mix together and sink along the 4°C isotherm' '(the thermal bar ) . The result 
is a four-cell circulation pattern with sinking motions at the thermal bars 
and upwellings in the middle and on both edges of the lake. Zonal flows are 
mostly restricted to the upper and lower Ekman layers^ The meridional cur- 
rents indicate loosely sheared flow along the thermal bar with the center of 
the lake as the nodal point. 

Some field measurements for currents and eddy viscosities were conducted 
during I967 and I968 in Lake Michigan. Comparisons have been made between 

xi 



these theoretical predictions and the previously published field data. Quali- 
tatively^ the thermal current structure and temperature distribution agree 
very well vith other observations as well as with our own field studies. This 
leads one to postulate that the thermal body force does play an important role 
in generating the mean lake circulation pattern. 
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I. INTRODUCTION AND LITERATURE SURVEY 



Introduction 



The St* Lawrence Great Lakes constitute the largest body of fresh water 
in the world and have long been considered as the inland "oceans" in the North 
American Continent. In spite of various differences between the properties of 
oceans and those of freshwater lakes^ general physical phenomena in the two 
remain quite similar. The special characteristics of a freshwater Great Lake^ 
such as negligible salinity and tides^ and the more tractable boundary condi- 
tions make the lake an ideal huge hydrodynamlc laboratory. In recent years 
great emphasis has been given to the value of the lakes as testing basins for 
oceanographic studies and to understanding all aspects of the limnological 
dynamics of the lakes. Figure 1 shows the Great Lakes and their mean surface 
currents during warmer months (after Millar^ 1952). 

The climate in the Great Lakes region is tjrpically humid continental with 
a severe winter and a warm summer. The continental climate^ modified by the 
stabilizing effect of the body of fresh water^ gives this region more rapidly 
changing and complex weather patterns than those of the maritime region. Since 
the Great Lakes lie in the interior of the North American continent between 
source regions of contrasting polar and tropical air masses^ this region is at 
the junction of the paths of storms from several areas of cyclogenesis in the 
western portion of the continent. Cyclonic storms produced along the polar 
front usually move toward the Great Lakes under the influence of the general 
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western circulation. The riximerous shifts in wind and changes in weather are 
due mostly to the many moving cyclones and anticyclones traversing the area. 
Although westerly winds prevail^ winds from any direction are likely to be en- 
countered, shifting in relation to the particular position with respect to the 
moving pressure centers (U, S. Dept. of Commerce, 1959) • 

The general annual thermal cycle in the lake has been discussed by Hutchin- 
son (1957)* For a relatively deep large lake in the temperate continental 
climate region, such as any one of the Great Lakes, the lake temperature falls 
below k°C during the winter and rises far above k°C during the summer. The 
lake is homogeneous in the winter and stratified in the summer. The quasi- 
homogeneous period in I^ake Michigan is from November through June (Ayers, 1962) . 
Fall turn-over continues almost all winter. As the season progresses the lake 
gradually becomes stratified with a relatively small vertical temperature gra- 
dient in the upper layer, the epilimnipn, and another relatively larger vertical 
temperature decrease in the lower cold layer, the hypolimnion. Between these 
two layers there is a very thin layer of sharp vertical temperature gradient 
called the thermocline. In the autumn the epilimnion thickens due to convective 
mixing induced by surface cooling, the thermocline becomes weaker and deeper 
and finally disappears, and the lake again becomes homogeneous* If the lake 
freezes in winter, inverse stratification may be set up below the ice. However, 
for the lakes in this region, the surface is only partially ice-covered during 
the winter. Therefore most of the main body of the lake cools progressively 
during winter and reaches temperatures below ^°C. The thermal cycle does not 
start until the arrival of the warm-up season. 



Since soil has a much lower specific heat and thermal conductivity than 
water^ the land surface is more sensitive to the ;rate of change of heating or 
cooling the atmosphere. Consequently^ a temperature gradient^ either lakeward 
or shoreward, always exists along the perimeter of the lake during warming-up 
or cooling- off seasons. Figure 2 shows the surface temperature pattern in 
July 1962 (McFadden and Ragotzkie^ 19^J>) * Figure 5 shows some typical surface 
temperature transects in Lake Michigan during warming-up and cooling-off periods. 
In the late summer and through autximn, there is not much horizontal temperature 
difference between land a.nd water and a strong thermocline always exists in the 
lake. In general, the lake is homogeneous, or nearly so, in winter, spring, and 
early summer. 

A unique density characteristic of water makes aquatic life in it possible 
in winter as well as in other seaspns. The density of fresh water is a maximum 
at k°C, and decreases as the temperature deviates, either positively or nega- 
tively, away from the temperature of maximum density. The variation of the 
density of water in relation to the temperature substantially affects the nature 
of heat distribution and develops circulation phenomena in the body of water. 

The energy inputs into the lake that are responsible for mixing and move- 
ment of its water, directly and indirectly, are the meteorological forces, the 
thermal forces, and the tidal forces. The last factor is negligible in the 
Great Lakes (Hough, 1958). Other forces which do not produce currents (but 
affect them) are the Coriolisj f rictional, turbulent, and centrifugal forces. 

Of the two main current generating forces, most scientists, in predicting 
lake currents, have considered the effects due to the wind stresses and neglected 
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the thermal forces. So far^ in my opinion^ no satisfactory result from trying 
to relate the general current circulation pattern of the lake to the wind 
stresses has ever been achieved. On the other hand^ though the temperature 
field in the lake is observed to be much more stable than the wind and far more 
detailed measurements of temperature data are available^ few people have tried 
to investigate the relation of the thermally induced forces and the lake cur- 
rents. 

Given the temperature distribution in the lake, are the thermally induced 
forces really so small that they can be neglected in modeling the lake circu- 
lation, especially under a constantly shifting and changing wind field? 

It is the purpose of the present study to model the flow in Lake Michigan, 
considering the geostrophic thermal gradient force as the prime source for the 
mean lake current generation, and to investigate the thermal circulation pattern 
of the lake that results. 

A review of the temperature distributions for all the stages of the annual 
cycle in the Great Lakes shows that the basic thermal phenomena in the lake can 
be approximately grouped into seven types: 

1. In the middle and late winter^ between January and early March, the lake is 
below 4°C everywhere. Due to differential cooling from both sides of the 
lake^ freezing starts from both shores and slowly progresses lakeward. The 
surface temperature of the partially ice-covered lake may be expressed as a 
cosine function with a maximum temperature, less than k°C, in the middle 
portion of the lake. This is the typical case for winter cooling. This 
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condition usually lasts for a quite long period of time until the lake is 
almost free from ice in the spring, 

2. In late March or early Aprils the warming-up season starts. During this 
period the whole lake is still below U^C and some small portions of the 
lake may still be covered with floating ice. Due to the quick response of 
the land to the warming of the atmosphere^ and the almost as rapid response 
of the shallow edges of the lake^ the surface temperature of the lake can 
be approximated as a small amplitude concave cosine function with maxima of 
about 4*^0 at the two edges of the lake. This is the period of spring- heating 

5. As the season progresses the temperature of the water near the shore is 
heated up to above ^°C while the middle part of the lake may well be at 
less than k°C. Then^ the special unique characteristics of water will in- 
duce a near-shore mixing process which has been observed in all large and 
deep lakes in temperate regions during the warming-up and cooling- off 
periods. 

During the warming-up period, the temperature of the inshore edges of 
the lake is all above k°Cf a.nd a direct temperature stratification exists 
in the edge waters, but the main body of water in the middle portion of the 
lake is less than ^''C and essentially isothermal. 

During the cooling- off period similar, but opposite, temperature con- 
ditions exist. The boundary band between these two temperature stratified 
regions, one isothermal and the other stratified, is called the "thermal 
bar." The thermal bar phenomenon is obviously connected with combination 
and mixture of waters from different sources at temperatures above and below 



the temperature of maxinrum density. The mixing of warmer and colder waters 
causes a. sinking motion along the 4°C isotherm^ and a type of isolation of 
the flows on either side of the ^° isotherm is established. Waters of dif- 
ferent color, marked differences in phytoplankton populations and different 
flow patterns are observed along the ^"^ mixing zone (Stoermer, 1968). Sink- 
ing along the k° isotherm requires surface convergence. The convergence 
implication of the thermal bar requires an upwelling in the middle portion 
of the lake during this period of time. This upwelling provides water 
colder than the temperature of maximum density at the horizontal center of 
the lake until the disappearance of the thermal bar. The spring thermal 
bar period usually lasts for four to six weeks. The surface temperature 
profile can still be approximated as a concave cosine function, but the mini- 
mum temperature at the middle of the lake is below 4°C. 
kn After the disappearance of the thermal bar the lake surface is above ^°C. 
The horizontal temperature difference between the edges of the lake and its 
middle portion has much increased and ascends to a maximiim. A positive 
temperature gradient toward the shore always exists in this stage and the 
surface temperature profile of the lake can be approximated as a concave 
cosine function with its minimum near ^°C at the middle of the lake, A 
weak thermocline starts to grow in late June or early July. At the begin- 
ning, a weak, shallow, thermocline may appear during the day and disappear 
at night. In general, the whole lake Is homogeneous, or almost so, during 
this period of time. This is the summer heating period. 
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5. After July^ the lake is strongly stratified. A sharp vertical temperature 
gradient thermocline is present all across the lake. The almost homogeneous 
epilimnion is nearly isolated from the cold stable hypolimnion. This is the 
strong thermal stratification period. 

6. The lake starts cooling off in late October or November as indicated in 
Figure 5, In December^ the lake finally becomes quasi-homogeneous again. 
At this time^ the whole lake is still above k°C but the edges of the lake 
have temperatures lower than the middle. This is the autumcn- cooling period. 
The surface temperature profile ret\irns to cosine form. 

7. After the edges of the lake have been cooled to a temperature less than 
^°C^ an autumn thermal bar phenomenon is expected to appear. However, the 
horizontal temperature difference between the edges and the main body of 
water, which is then nearly ^°C, is much smaller than the spring warming-up 
period. Therefore, the cooling-off thermal bar phenomenon is much weaker, 
and may even hardly be recognized. 

After the cold weather becomes severer, the whole lake will soon be less 
than ^°C but is still cooling. Then the temperature cycle starts again* 
In our study the following circulation patterns are obtained: 
Type -A Circulation : During the period of winter -cooling and of siimmer- 
heating, a cylconic circulation pattern of meridional velocity (north- south flow) 
is found in the lake with the middle of the lake as the nodal line. The zonal 
component (east-west flow), together with the vertical component of the velocity 
field, forms two cells of circulation in the vertical cross- section of the lake. 
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There is a general weak"^ broad sinking motion in the middle portion of the lake 
and upwellings near both shores. Two narrow Ekman layers having zonal flows on 
the free surface and at the rigid bottom connect the cell circulations. 

Type -B Circulation ; The circulation pattern d.uring the period of autumn- 
cooling and of spring- heating time Is found to be almost opposite to that of 
the type-A circulation. 

Type -C ■ Circulation : This i$ the thermal bar circulation pattern^ either 
of the spring warming-up period or of the autumn cooling- off period. The cir- 
culation pattern is more complicated^ as expected. The meridional velocities 
form two sets of weakly sheared flows along the thermal bars on each side of the 
lake. In the cross- section of the lake., the, zonal and the vertical velocity 
components form four cells^ two broad and large cells in the main body of water 
in the middle portion of the lake between two thermal bars^ two other relatively 
more intense cells between the thermal bars and the shores. The cells circulate 
clockwise and counterclockwise in. such a way that there is a general weak up- 
welling in the middle portion of the lake^ sinking motions along the thermal bar 
mixing zones^ and upwellings near the east and west shores. 

The current circulation pattern for the strong thermal stratification 
period is not predicted in a quasi-homogeneous fluid model study. 

To verify the importance of the thermal forces^ it is convenient to compare 
the real field data with the mathematical model prediction in the corresponding 
periods of time. The temperature distribution and the velocity fields from the 
theoretical solutions qualitatively agree very well with field observations. 
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This leads one to postulate that the thermal force does play an important role 
in generating the mean lake circulation patterns. 

Previous ObgejnratJQnal ' Studies 

Lake Michigan is the third largest of the Great Lakes^ the only one totally 
within the United States. It is also one of the most surveyed large lakes in 
the world. For convenience, we are going to use the physical setting of Lake 
Michigan to investigate the typical thermally induced current structure in the 
Great Lakes. 

The first study of the surface current in the Great Lakes was conducted by 
Harrington in the period I892 to l89i^- (1895). Figure 1 shows the general cur- 
rent patterns during warmer months found by Harrington and drawn by Millar 
(1952). Harrington's results from drift bottles show that the main surface 
current in Lake Michigan circulates in a cyclonic direction around the lake 
together with two separate cells, one in the southern basin, the other in the 
northern basin. Deason (1952) studied Lake Michigan currents using the same 
method and found a "swirling" circulation whose direction he does not give. 
Church (I9^2a, 19^2b) made a comprehensive study of the annual thermal struc- 
tures in many cross- sections of Lake Michigan and was the first to try to relate 
the currents to the observed temperature distribution. His remark, stating that 
the denser water is in the middle of the lake d\H*ing most of the year, implies 
a cyclonic circulation in Lake Michigan. Ayers et al. (1958) made a series of 
synoptic surveys in the Great Lakes including four cruises in Lake Michigan, 
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studied the currents "by using drift bottles, and computed the currents in Lake 
Michigan using the modified dynamic height method (Ayers/ 1957). Their June 
currents showed a cyclonic circulation but the August current had a different 
pattern. Qualitatively, Ayers^ drift bottle data agree well with Harrington *s 
study. In 1962 and 196^, in order to understand the water pollution problem in 
Lake Michigan, the Federal Water Pollution Control Administration (FWPCA) under- 
took a lakewide long period detailed survey of the temperature and current by 
deploying recording instruments all over the lake at thirty to forty stations « 
At each station. Woods Hole type temperature recorder and Richardson current 
meters (Richardson, I962) were set as pairs at fixed depths: 10, 15, 22, and 
30 m, and at each succeeding 30 m level. A wind anemometer was installed on the 
buoys of some of the stations. Drogue studies for water movement were also 
made during that period. "This is the first occasion in which physical events 
and processes in so large a natural water body, fresh or marine, have been 
monitored for an interval of nearly two years with so close a network of record- 
ers" as Mortimer remarked in the foreword of the report (U. S. Dept. of Inter- 
ior, 1967). Unfortunately, most of these long series of recorded data have not 
yet been analyzed and published. From the analyzed data of some stations, the 
report indicates that, generally speaking, surface currents are weak and south- 
bound near Lake Michigan's western shore, but narrower, stronger and north- 
bound near its eastern shore. No positive correlation between the wind above 
and the surface current at 10 m depth could be found. No permanent type of 
circulation pattern exists in the lake, but four basic patterns have been ob- 
served« One winter pattern shows the currents circulate counterclockwise with 
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gyres in the northern and southern basins ♦ The principal flows are northward 
along the east shore and southward along the west shore. This pattern is most 
often observed between January and April.. The other winter pattern shows re- 
verse circulation and is usually found after November. The summer patterns 
have more complicated pictures. The prevailing pattern shows the existence of 
a cyclonic gyre in the southern basin, but inshore currents on both shores are 
northward, in spite of the fact that the northern flow along the west shore is 
much weaker. In the northern basin, the dominant flow is southward in the cen- 
ter of the lake. The other summer pattern still has the cyclonic gyre in the 
southern basin but the inshore currents along both the east and the west shores 
are moving southward. The subsurface current has no permanent circulation pat- 
tern. Though the cyclonic pattern prevails, there are many instances in which 
the pattern is reversed. The day-to-day movements suggest that there is as 
much variability at 2^0 m as to the 10 m level. During the period of strong 
thermal stratification, it is found that the currents above and below the thermo- 
cline move in opposite directions. The studies also showed that currents exist 
through the entire vertical column of water and the upper layer is mostly wind- 
driven. Mortimer (1965) postulated that internal oscillations in the form of 
Kelvin waves near the shores and Poincare waves with periods less than inertial 
period (about I7.5 hr in mid-lake) in the middle lake are responsible for the 
current circulation pattern in Lake Michigan during the period of strong stra- 
tification. Csanady^s (I967, 1968a, 1968b) studies of a simplified two-layer> 
constant depth, "model Great Lakes" with wind-driven currents, have suggested 
the existence of the following phenomena during the summer stratification: a 
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"coastal Jet^" large thermocline movements near the shores^ surface and inter- 
nal oscillations, and rotary current in the center portion of the lake* Verber 
{I96hf 1965) has found the presence of inertial rotary currents and internal 
waves for prolonged periods during the stratified period in Lake Michigan, 
Noble (1966), using FWPCA Great Lakes- Illinois River Basins project data (sta- 
tion 8), attempted to synthesize the observed current structure with the Ekman 
type wind-driven model (Ekman, I905) and was not successful. Alternatively, he 
suggested a geostrophic eddy model with the existence of contrarotating eddies 
in the lake. The persistence of the surface temperature profiles indicates 
stability of the geostrophic circulation (Noble, I967) . In studying the ther- 
mal structure in the Great Lakes, Rodgers (1965) has described the thermal bar 
feature in detail. Field data during the spring from all the Great Lakes con- 
firm Tikhomirov^s (1965) suggestion that the thermal bar as was observed in 
Lake Ladoga is a common phenomenon to all large dimictic lakes in the temperate 
zone. 

Literature Survey 

In geophysical fluid dynamics the Coriolis force plays an important role 
and must be taken into consideration. In Lake Michigan, except during the 
strong thermal stratification period, a horizontal temperature gradient always 
exists, especially during the warming-up and cooling-off seasons as shown in 
Figures 2 and 5. In the literature of differentially heated rotating flow, 
Hadley (1755) was the first to notice that the longitudinal atmospheric circu- 
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lation is the result of the body force induced by the latitudinal variation of 
the heating balanced by the Coriolis acceleration. In an attempt to understand 
the general circulation of the atmosphere and the origin of the earth's magnetic 
fields^ a series of experiments on rotating differentially heated fluid, con- 
tained either in a "dish pan" or an annular container, were conducted by Fultz 
(1955)/ Hide (1955^ 1958)^ Fowlis (I965), and Fowlis and Hide (I965). If the 
viscosity of the fluid is relatively small and the angular rotation is moderately 
high, the experimental results show that four distinct regimes of hydrodynamical 
flow may appear. The principal properties of the flow depend mainly on certain 
nondimensional parameters, such as the thermal Rossby number and the Taylor 
number. Figure k, reproduced from Barcilon (1964), shows the transition curve 
obtained from Fowlis* experimental results (I965) on the (h) , Ta plane in which 
all flow regimes are observed. The parameter 

® " n^(b - a)^ 
is the thermal Rossby number, and 

Ta = ^ B. ^ 

V d 

is the Taylor number, where b and a are the outer radius and the inner radius 
of the annular container j, respectively, d is the depth of the fluid in the 
annulus^ a is the thermal expansion coefficient^ v is the kinematic viscosity, 
AT is the impressed temperature difference, g is the gravitational constant, 
and Q is the rotation rate. As shown in the figure ^^ there exists a value of 
the Taylor number (Ta)^ say Ta*^^ below which axisymmetric flow arises irrespec- 
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tive of the value of the thermal Rossby number ((H)). The upper or lower axi- 
symmetric regions also occiir when Ta is greater than Ta^ and @ is sufficiently 
high, greater than a particular value (H) i, or (9) is sufficiently low, smaller 
than another particular values (3) 2* The nonsymmetric regimes appear only when 
Ta is greater than Ta* and (S) has the value between @^ and @2- ^^^ 
nonsymmetrical flow pattern may be in regular wavy form with jet- streams, or, 
after some transitional vacillation phenomena, the flow may turn into irregular 
complicated fluctuations. Theoretical stability analyses on this problem have 
been carried out using different approaches, by Eady (19^1-9)^ Davies (195^, 
1959), Kuo (1957), Brindley (i960), and Barcilon (196^). Robinson (1959), and 
later Hunter (1967), investigated the lower symmetrical regime for fluid con- 
tained inside a rotating annulus of square cross- section subjected to horizontal 
differential heating. An entirely geostrophic thermal wind is obtained over 
the maiin body of the fluid in the limit of small thermal Rossby numbers. In the 
cross-* section of the annulus, the circulation is much weaker and exists only 
within boundary layers circumscribing the rigid walls. A large cell is found 
to occupy most of the cross- section and two additional small cells are confined 
to the side-wall boundary layers. If the upper surface is free, the vertical 
side boundary layers have a double structure with one boundary layer inside 
another (Hunter, I967) . 

Hide (l96Jt) has studied the viscous boundary layer at the free surface of 
a rotating baroclinic fluid and found that the ageostrophic flow in a free 
surface boundary layer is much weaker, by a factor of the square root of the 
Ekman number, than the corresponding flow in the layer at the rigid surface. 
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In the thermal structure problem, in spite of the fact that the free surface 
Ekman layer is weak^^ vertical heat advection induced by boundary layer suction 
at the free surface can dominate thermal conduction if the product of the 
Prandtl number and the thermal Rossby number is not smaller than unity. (The 
corresponding result for a rigid surface is PR « e^/^.) 



II. A BAROCLINIC DYNAMIC MODEL OF LAKE MICHIGAN 

General Physical Description , of Lake Michigan ^ " 

Lake Michigan has a surface area of 58^016 km^ and a total volume of W78 
km^. The average depth of the lake, including the shallow Green Bay is 84 m, 
with the maximum depth of 281 m in the northern basin^ The length of Lake 
Michigan, over i^-94 km, is much greater than the mean width of 120 km. Except 
for a few irregular points, most of the 2673 km shorelines are straight and 
smooth, and nearly in north-south orientation. An east-west cross-section of 
the central part of Lake Michigan from Milwaukee, Wisconsin, to Muskegon, 
Michigan, is shown in Figure 5. The average beach slope around the lake is 
about 1:100 (Davis and McGeary, I965 ) • 

Lake Michigan and Lake Huron are, hydraulically speaking, only one lake. 
The connection between these two lakes is the broad and deep Strait of Mackinac. 
Therefore, their surfaces stand at the same level of 176,5 m above sea level* 
With the exception of some water diverted at the Chicago Canal into the Illinois 
River, the entire discharge of the Michigan basin goes to Lake Huron through 
this strait. Water budget estimates show a net outflow at the Strait of 
Mackinac of lljii m^ per second (U.S. Dept. of Interior, 1967)- Using this av- 
erage value, the total annual out-flow is calculated to be much less than \io 
of the total Lake Michigan volume. Therefore, the hydrological effect is not 
very important . 
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The Coriolis force per unit volume flow is twice the normal component of 
the earth's rotation vector at a particular location. Therefore^ it varies 
with latitude. In Lake Michigan its magnitude ranges from O.976 x lO"""^ rad/sec 
in the southern end to I.05 x 10""^ rad/sec in the northern end. The variation 
from the south end to the north end is small. Hence^ the beta-effect in Lake 
Michigan is not important. 

It is an obvious and well observed fact that wind produces waves ^ induces 
currents, and results in turbulent mixing. Ekman's theory (Sverdrup et al., 
19^2) on the wind-driven current in the open sea has been verified and doc- 
umented. Shallow water drogue studies (about I.5 m deep) in Lake Michigan 
show the surface currents are mostly wind-driven with a response time constant 
ranging from one hour to several hours after a wind shift. The direction of 
surface currents differ from that of the mean prevailing wind by from 20° - 80® 
(U.S. Dept, of Interior, I967)* However, various ways of weighting past wind 
averages used in the FWPCA project in trying to correlate the wind and the cur- 
rent water data 10 m below in Lake Michigan, show no definite correlation (U.S. 
Dept. of Interior, I967). 

In the Great Lakes area the wind system is modified by the continuous 
presence of a mesoscale high pressure in the spring-summer period and a low 
pressure in the fall-winter period due to the meteorological lake effect. The 
wind observed on the lake is quite different from that of nearby land stations. 
Hunt (1959) indicates that a' wind differential exists between the land and the 
lake. In Lake Michigan, the winds over the lake may be 96^0 greater than those 
over the city of Chicago at certain times of the year (U.S. Dept. of Interior, 



25 

1967). As stated in the Introduction, the wind in the Great Lakes region has 
numerous changes both in magnitude and in direction and is seldom constant over 
a long period pf time^ The apparent winds may well be felt differently at 
different locations in the lake even if there is no change in the wind system. 
Though wind is responsible for inducing the surface current^ the long response 
time of the currents at deeper layers may well average out the effects of wind 
stresses over a long period of time, for rapidly changing winds. Therefore, 
as far as the general mean lake current circulation is concerned, it is possi- 
ble that the wind-driven current may npt be the most important. This is espec- 
ially obvious during the lake warm-up season. In the spring, when the lake 
water is still relatively much colder than the already warmed adjacent land^ 
the air mass above the lake is frequently in a stable condition. Wiresonde 
cross-sections in Lake Michigan have demonstrated the intense spring time in- 
version prevalent before the warming of the water surface (Bellaire, I965 ) . 
At the Wisconsin shore, it Uas be^n observed that water remains smooth even 
beneath offshore winds of greater than 10 knots (at I6 m) (Strong and Bellaire, 
1965 )♦ Therefore, during this period the meteorological force is even less 
important. However, in general, turbulent mixing in the open lake is mostly 
due to wind stresses. Wind is also responsible for set-ups, seiches, and sur- 
face waves. These effects may have a secondary influence on the water movement 
and mixing process in the lake. 

The temperature structure in Lake NtLchigan is in a continual quasi-steady 
state. The temperature is, generally speaking, a slowly varying function of - 
time. During the warm-up season the sun and river runoff are putting energy 
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into the lake longer and more consistently than does wind stress. There are 
some changes in magnitude but solar energy heats the lake during the whole sea- 
son. And it is similar, with diminishing trendy for the cooling-off season. 
During the lake warming-up and cooling-off season, temperature gradients, 
either horizontal or vertical, show persistence, which, in turn, suggests that 
the thermal energy input may dominate the lake circulation dynamics at this 
time. This is evident from the persistence in the thermal features in the 
transects as shown in Figure 3 (Noble,. I966) . 

The frictional forces are important near all boundaries. In a flow field 
of a turbulent nature, as in Lake Michigan, the frictional forces are expressed 
as the space rate of change of the products of eddy viscosity and Reynolds 
stresses, in analogy to laminar flow. The average horizontal eddy viscosity 
in Lake Michigan is the order of 10^-10^ cm^/sec and vertical eddy viscosity 
(l-lO^ cm^/sec) is much smaller. Eddy diffusivity may be of the same order as 
eddy viscosity but it is a bit smaller in general. Eddy viscosity and eddy 
diffusivity of the fluid depend not only on temperature (such as the molecular 
viscosity), but change from place to place and time to time, depending on the 
local character and dimension of the flow. The overall average values of the 
order of magnitude of eddy viscosity and eddy diffusivity can be crudely treated 
as "constant" for application . to the present problem (see Chapter V). 

In general, though a weak thermocline may exist during the early summer^ 
heating and late autumn ■ cooling period, the lake is qua si -homogeneous during 
the warming-up and the cooling-off seasons. Of course, in the late summer and 
early autumn season^ the lake is strongly stratified. 
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Modeling Considerations 

During the warming-up and cooling-off seasons in Lake Michigan, a sharp 
horizontal temperature gradient always exists along the perimeter of the lake. 
Given the variability of the wind fieldj, it is reasonable to assume that the 
time-averaged current is due to thermal forces. Equilibrium between the thermal 
body force and the Coriolis acceleration produces a horizontal velocity per- 
pendicular to the direction of the existing thermal gradient. This is the geo- 
strophic thermal wind relationship, which is considered to be of dominant im- 
portance. 

The length of Lake Michigan is much longer than the average width. The 
north end of the lake is connected through 9 wide and deep strait to an even 
bigger lake* In the Lake Michigan model, to the lowest order of approximation, 
the flow field, away from the southern end, can be considered as independent 
of latitude* 

The primary current generating force under consideration, the thermal body 
force, is a slowly varying function of time. In studying the thermal current 
structure of the lake for different seasons individually, all variables can be 
considered as time independent without^ losing good qualitative approximation* 

The tremendous difference between the horizontal eddy viscosity and the 
vertical eddy viscosity makes the choosing of different values for the viscos- 
ity a necessity* But, in general, the vertical velocity gradient is much 
greater than the horizontal. Therefore, the effect of mixing due to the ver- 
tical transport may well be of the same order as the lateral mixing* In con- 
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sidering the horizontal Reynolds stress of the same order as the vertical 
Reynolds stress, the ratio of the vertical to the horizontal eddy viscosity 
is in the same ratio as the depth to the horizontal width of the lake, which, 
incidently, agrees well with our field data. 

In our preliminary model, we will not consider the strong thermal strat- 
ification period of the year. 

The Boussinesq approximation is used in the formulation. Since the largest 
horizontal dimension of the lake is much smaller than the radius of the earth, 
curvature effects can be neglected. The centrifugal acceleration is neglected 
as compared with the gravitational acceleration. In the energy equation, the 
dissipation term is considered as small and therefore is neglected. 

Mathematical Formulation of the General Governing Equations. 

Together with the above restrictions, let us assume that the fluid is con- 
tained in a long symmetrical "trough," trapezoidal in cross-section. It is 
subjected to a surface temperature with a maximum or minimum at the center. 
The lake center is taken as the center of the Cartesian coordinates (Figure 6). 
The lake is rotating with angular velocity (component of the earth's rotation 
normal to the lake surface) about its vertical axis, which is anti-parallel to 
the gravitational force. The horizontal surface of the lake is considered free 
from wind stresses. The dominant flow is assumed to be meridional, thus the 
basic geostrophic thermal gradient relation is satisfied. All specifications 
of nomenclatures are summarized in Table 1» 
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-l/2+X1(l-z) 



Figure 6, Cross-section of a long, symmetrical, 
trapezoidal Lake Michigan model. 



The equations governing the steady motion of an incompressible, viscous^ 
heat-conducting, rotating fluid are: 



q*-Vq' + 2^ X q^ = - 



— V^p' 

^o 






'e ax' 



Yhz' 



V'-q' = 0, 



q'-V'T- = (K„-^ + K,,|rr^) T- 



H ax' 



vaz' 



(1) 



(2) 



(5) 



which are, respectively, the Navier-Stokes equations, the continuity equation, 
and the energy equation. Note that the above dependent variables are inde- 
pendent of y. V is a two-dimensional gradient operator. All primes are di- 
mensional quantities. 
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TABLE 1 
Nomenclature and Nondimensional Parameters 



Symbol Definition 



X Horizontal east-west coordinate^ zonal direction 

y Symmetrical north-south coordinate^ meridional direction 

z Vertical coordinate^ negative is in the direction of gravity 

q Velocity vector with components (u^ v^ w) 

p Deviation of the pressure from hydrostatic pressure 

T Temperature 

(N.B. The above symbols primed have dimensions^ unprimed are dimen- 
sionless . ) 

T Reference temperature, ^°C 

A Constant for density approximation 

(i^j^fi,) Unit vector in the (x^y^z) directions 

g Acceleration due to gravity 

a Coefficient of thermal expansion; always positive except below ^°C 

V Kinematic eddy viscosity, subscript H denotes horizontal, V, vertical 

K Thermometric eddy diffusivity. Subscript H denotes horizontal, V, 
vertical, 

£1 Angular velocity of the vertical component of earth's rotation in 

Lake Michigan 

L Width of the lake, reference horizontal length dimension 

D Depth of the lake, reference vertical length dimension 

AT Total horizontal temperature difference between the edge temperature 
and the reference temperature T 

■^ Stream function for the zonal and the vertical velocities 
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TABLE 1 (Concluded) 



Symbol Definition 



p Fluid density at temperature T 
7 Depth to width ratio 

K Tangent of the co-slope angle with the vertical^ \ = tan 9 (see 

Figure 6) 

V Two-dimensional gradient operator 

^ — + ic — 



V^ Two-dimensional Laplace operator 
^ 7 Sz 



^ A fourth order linear operator 

SI ^ . ^ ^ . ^ ; ^ 1 a^ 



2 



V Dimension of velocity 

V = agDAT/2nL; or V = AgDAT sin(AAT) 

2 nL 

R Nondimensional thermal Rossby number 
OgDAT ^ AgDAT sin (A AT) 



Nondimensional Ekman number 

^ 2J^LD 
Eddy Prandtl number 

K 
V 
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Boundary Conditions and Seasonal Variations of the Surface 
Temperature Distributions 



VELOCITY FIELD 

The boundary conditions for the velocity field are the no slip conditions 
at rigid surfaces, the bottom and the slanted side boundaries, and the condi- 
tion of zero stress and zero normal component of velocity at the free surface. 
The latter is assumed to be plane and horizontal. Hence, we must require that: 

SEASONAL VARIATION OF THE SURFACE TEMPERATURE DISTRIBUTION 

As discussed in the introduction, during the winter cooling period, when 
the whole lake is below k'^C and both edges of the lake are colder than the 
middle, the surface temperature profile is in the form of a cosine function 
with its maximum, slightly less than 4°C, at the middle. The current circula- 
tion induced in our model during this period has the same pattern as in the 
summer heating period. During the summer heating period the surface tempera- 
ture in the lake is above k°C and the edges have much higher temperatures than 
the middle. The fact that the same circulation patterns are induced by these 
two imposed temperature forms is due to the special characteristics of fresh- 
water density and of the change in direction of the thermal gradient during 
the heating and the cooling periods. This type. of circulation pattern is called 
the type-A circulation. In the mathematical demonstration of the type-A cir- 
culation pattern, we are going to use a concave cosine function as the free 
surface temperature condition. That is: 
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T^ := ~ (1 - COS ^Y^) AT + T at z' = D. (5) 

The type-B circulation pattern is expected to exist during the period of autumn 
cooling when the whole lake is still above ^**C with only two edges at a lower 
temperature near i^-*'C, and during the period of spring-heating when the whole 
lake is all below U°C with edges in a bit warmer state at nearly J^°C, The 
type-B circulation is the exact opposite of the type-A circulation. 

The type-C circulation pattern is the current circulation during thermal 
bar periods. The surface temperature of the spring thermal bar period is 
characterized by the existence of a sharp temperature gradient between the much 
warmer boundary water and the colder main body of water in the middle. In the 
middle portion of the lake, the temperature may be well below i+°C but the two 
edges have reached a temperature far above U°C. Therefore, letting the thermal 
bar be near the temperature of maximum density, we require 

T' = 7 (1 - 2 cos ^-) AT + T at z = D (6) 

5 L o 

where AT is the temperature difference between the maximum water temperature 

near the shores and the thermal bar temperature. The autumn thermal period 

has almost the opposite temperature profile except AT is smaller and negative. 

The lake temperature near the bottom is very stable all year round with 

little variation from the temperature of maximum density. Therefore, we assume 

the bottom temperature to be constant at nearly 4°C in all cases. As to the 

slant side boundaries, we require the insulation conditions to be satisfied at 

X = ± L/2, which permits some conduction between the earth of the slant shores 

and the edges of the lake water. This makes the side boundary conditions even 
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more similar to the real situation in the lake. 

The differential equations (l)^ (2)^ and (5) together with all the bound- 
ary conditions form a high order, coupled, nonlinear system. Exact analytic 
solutions are very difficult to obtain. Therefore, some approximate methods 
capable of making the system more tractible are preferred. In the next chap- 
ter, the general thermal currents, including type-A and type-B circulations 
will be investigated mathematically in great detail. 



III. GENERAL THERMAL CURRENTS 

Mathematical Formulation and Methods of Appr oximat ion 

THE EQUATION OF STATE 

Except during the thermal bar time, the temperature in the lake is either 
all above k'^C, as in the summer -heating and the autumn-cooling periods, or all 
below k°C as in the spring -heating and the winter-cooling periods. Then, using 
the Boussinesq approximation, the equation of state is a linear relation be- 
tween the density and the temperature as 

P* = Po^l - C^(T' - T^)] . (7) 

where a, the coefficient of thermal expansion, is treated qualitatively as a 
constant in each individual period, and T^, the reference temperature, is k^'C* 
Hence, the thermal buoyant force in the momentum equation is: 



^ = - a(T^ - To) (8) 

Po 



This approximation is valid for all the type-A and type-B circulations, which 
are called general thermal currents to distinguish them from the thermal bar 
circulations. 

NONDIMENSIONAL GOVERNING EQUATIONS 

It is convenient to introduce nondimensional variables and parameters. 
The following unprimed nondimensional variables are defined: 
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5^ 



T' - T 
Lw' _, n' 

w = — — and p 



(9) 



DV ^ p agDAT ' 



where V is the reference dimensional velocity and all other variables have been 
specified in Table 1» The vertical velocity is nondimensionized in this form 
in order to satisfy the continuity requirement « 

In our basic assumption^ we consider that the dominant mean flow field is 
in geostrophic thermal gradient equilibrium* Then the basic thermal "wind re- 
lationship should be satisfied^ that is 



-i = <*i- (-) 



This gives V^ in dimensional form^ as; 



V = ^2AI • (11) 

2QL 



where D is the depth of the fluid* Since AT is the horizontal temperature dif- 
ference^ and the baroclinic circulation is induced by the horizontal density 
difference^ the reference velocity V is the anticipated order of magnitude of 
the meridional motion. 

In considering the mixing effect of vertical transport to have the same 
order as the lateral mixing effect^ the horizontal eddy viscosity to the ver- 
tical eddy viscosity has the same proportion as the reference length L to the 
reference depth Da That is^ 
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^H " D \ • ^^2) 



Equations (9), (H)^ and (12) are substituted into (l), (2), and (3). The 
nondimensional continuity equation is still the same form as equation (3) which 
permits us to define a stream function \i; by expressing the nondimensional veloc- 
ity as: 



q = jxVt(x,z) + vj . (13) 



The vorticity equation is obtained by eliminating the pressure by cross-differ- 
entiation from the X and z momentum equations. After using equation (I3), it 
becomes^ 

^^ Sz ax ^ |_az axdz^ ax az^ ^ ^ax ax^az az ax^-'J • ^^^> 

The y momentum equation^ noting the y-independence of the model, becomes^ 



^-t- «(||-||) • (^5) 



The eddy diffusivity is assumed to be of the same order as eddy viscosity^ 
then equation (3) leads to: 

^^ - -<! I - 1 1) • (^«) 

All nondimensional parameters and operators are specified in Table 1. Notice 
that: 
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QgDAT 
(2aL)^ 



is the nondimensional thermal Rossby number similar to the parameter(H)in Fig- 
ure k; . 



2fiLD 



is the nondimensional Ekman number which is similar to the reciprocal of the 
square root of the Taylor's number, and 

a = — y 

V 

is the eddy Prandtl number "which is of the order of unity« 

METHODS LEADING TO THE APPROXIMATE SOLUTIONS 

We are interested in the range of small thermal Rossby numbers and small 
Ekman numbers « Since R « 1^ examining equations {ik) and (15)^ we find that 
the limiting case is that of geostrophic thermal gradient force domination in 
the flow field _, "which is the basic assumption of the modeling* This leads to 
the method of regular perturbation by expanding all dependent variables in 
powers of R. 

The expansion of temperature in a power series in R needs Justification* 
From the energy equation (l6)^ since o is of order unity^ R totally determines 
the relative importance of the heat conduction or the heat convection. Small 
R means heat convection is less important^ However^ the Ekman number, e, is 
an even smaller number. At first inspection of equation (l6), it seems that 
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the temperature gradient terms coupled >/ith momentum advections cannot be neg- 
lected even in the zeroth order of approximation because the combination of 
these three parameters (aR/e) is apparently not a small parameter. Neverthe- 
less^ due to the anistropy of the induced zonal and vertical velocity fields 
the series expansion in powers of R for the temperature is valid everywhere 
in the flow fields Because the free surface boundary conditions force the 
stream function to be sraall^ of order e throughout the field, this makes heat 
convection as represented by the right hand side of equation (l6) still of 
order R. At the boundaries, as we will see later, since all velocity compo- 
nents and all high order temperatures do have boundary layer characteristics, 
the contribution due to the velocity coupled heat convection is still small, 
even smaller than that in the interior region. Hence, the relative unimport- 
ance of the coupled heat convections is made clear. Actually, as will become 
obvious later, leading terms of every order of the R-expansion of temperature 
are of the order unity. The expansion of temperature in powers of R is at 
least apparently asymptotic. 

In the regular perturbation, we express all dependent variables in powers 
of R as 

^ = Jo h"^'"' (17) 

and substitute into equations (li|-), (I5), and (l6)» After equating each power 
of R, we have the zeroth R-order equations as: 

y (o) bv^°^ dT^°^ 
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,v2^(o) _ |i 



(o) 



hz 



W°^ = 



(19) 
(20) 



Notice that the zeroth R-order temperature is induced entirely by conduction. 
The first R-order equations are : 



,^,M,^"^ ^^"^ 



bz Sx 



hz BxSz^ Sx hz^ 



-jH 



^>^ ^V°l di^°^ W°'^ 



9k SzSx Sz Sx' 



(21) 



eV^- 



(i) at 



(1) 



az 



9z Sx Sx Sz 



(22) 



5z 9x ^ Sz 



(23) 



The second Reorder equations are: 



oz ax 



at 



!!la!till 
axaz^ 



ax sz^ 



Sz SxSz^ Sx Sz^ 



+ 7^( 



at^°^ a-t^'^ 



dz Sx' 



(2^ 



a/!l aft^ , a£l a!i^ V^aft^. 

Sx SzSx^ Sz Sx^ " Sx 9zSx^ ^ 
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Sz 5z Sx Sx Sz dz 3k Sk Sz ' 

(25) 



,W°) = o(#|^-#|^.|*^fE^-f^|t^). (36, 

dz dx OK dz dz dx dx dz 



We are going to use second R-order equations for estimations of the approxima- 
tion error. With reference to the last chapter^ especially equations (l^) and 
(15)^ boundary conditions for the zeroth R-order equations are: 

(o) Si^^^ (o) (o) , , 

r = ^ = v^ ^ = , T^ ^ = at z = , (27) 

dz 



^(o) ^ |f|^ . |v^ = 0, t(°) = i (1-cos 2nx) at z = 1 , (28) 
oz oz 2 



(o) , |tl!l , ^(o) , „ 
dx 



on the slant side boundaries and 

^^°^ 1 

^- = at X = ±- . (29) 

For all high R-order equations^ the boundary conditions for the velocity field 
and the temperature on the slant boundaries are exactly the same as those of 
the zeroth R-order equations^ only temperatures on the free surface and the 
rigid bottom require: 

T^^^ = , i = 1,2,3^ ••• etc. , (30) 
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for higher i-th R-order temperature* 

In all the R-order equations^ the Ekman number e appears only in the 
differentiated highest order terms* Since e is a very small number connected 
with viscosity^ the highest order differentiated, terms may have contributions 
only in the boundaries where they are large enough to overcome the smallness 
of £• Therefore^ it becomes a singular perturbation problem* The method, of 
matched asymptotic expansions^ or generally^ the method of boundary layer anal- 
yses can be used (Van Dyke^ 196^4-^ Carrier 1955)* Following the procedures 
of the boundary layer additive type analysis^ let^ 

f = ^f + f J - 1,2,3,4 . (31; 

where f is any dependent variable; and the right superscript i denotes the i- 
th R-order in the series expansion of the variable f ; the left subscript j 
denotes the i-th R-order of the variable in the j-th region, j = 1,2,3,4; and 

the subscript o represents the overall part of the variable f, including the 

(i) 
interior region* Therefore, f has no boundary layer characteristics 

o 

.f \, j = 1,2,3,4, represents the boundary layer part of the variable in its 
J 

respective boundary region* If J = 1, that is the free surface boundary re- 
gion, d = 2, the rigid bottom region, j = 3 and j = 4 are respectively the 
eastern and the western slant boundary regions* Figure 7 shows the interior 
and all boundary regions of the Lake Michigan Model* For example, the zeroth 
R-order stream function^ t , is represented as 

(o) (o) (o) (o) (o) (o) 



in 



(o) 
The \|f is dominant in the interior region and also reaches all boundary re- 
o 

(o) 
gions without changing its order of magnitude, it is significant only on 

the free surface boundary and vanishes exponentially in the interior. Sim- 
ilarly, st y 3^ } 4^ Q'^re important only at their respective boundaries. 

(o) 
In the interior region, only the \|r part is important. 

ah 



(0) 




(2/ 




(0) indicates the interior region 

(1) indicates the surface region 

(2) indicates the bottom region 

(3) indicates the eastern boundary region 

(4) indicates the western boundary region 

The shaded areas are comer regions 

Figure 7. Interior and boundary- layer regions of the Lake Michigan Model. 

BOUNDARY LAYER SCALING 

Equations (18) and (19) are coupled and show boundary-layer characteris- 
tics. In boundary-layer analysis the boundary layer thickness can be esti- 
mated by equating the highest-order differentiated terms connected with the 
small viscous parameter to the other terms with no connection to the viscous 
parameter in the same equation. Therefore^ we shall evaluate the order of 
the boundary layer thickness as a function of the small Ekman number e* 

Equations (18) and (19) together with equation (20) imply tha^ 

,=^^,(°) , ||^ . , (53) 



k2 



where 



^^ = t^C^M=.*x)^Mi.f)^.ig,; 



is the sixth order linear operator* 

On the free surface and the rigid bottom we may scale the normal coordi- 
nates near the boundaries as 

C = el^^(z - b) , (33) 

and with no change in the tangential coordinates^ where b = on the bottom 

and b = 1 on the free surface* 

At the boundary^ the frictional forc^ is important and in equation (32) 

the products of the viscous parameter^ e^ and the leading terms in the sixth 

order linear operator may have the same order of magnitude as the other term 

with no viscous interaction. Substitute equation (55) into equation (52) and 

equate the leading viscous term to be the same order in 8 as the non-viscous 

term* This leads to 61 = 1/2 and e^ = e/y* Then^ with the conditions that 

the boundary layer part of the zeroth R-order stream functions^ .\|f , J = 1 

J 

or 2y is of only local importance and vanishes away from its respective bound- 
ary^ the equation (52) can be integrated twice and yields 

where j = 1 on the free surface and. j = 2 on the rigid bottom. 

The eastern and western slant side boundaries are represented^ in non- 
dimensional form by 



h-5 



X = ± 



[l - ^.7(l - z)] 



where \ is the co-slope of the slant boundaries^ that is^ \ = tan 9. As shown 
in Figure 6^ the co- slope angle^ 0, is the physical topographic average angle 
of the beaches inclined with the vertical line. Let 



i - s^^ 



l^x + I - \7(1 - z)J , 



(35) 

T] = Z 

be the local coordinates normal and tangential to the slant boundaries. Sub- 
stituting equation (35) into equation (52), with all the highest order normally 
differentiated terms combined in making the same contribution as the other term 
in the equation, will lead to 62 = I/2 and e^ = eK, where 

r -| '■/^ 



is a constant. Then, equation (52) becomes exactly the same form as equation 
(54.) except j = 5 for the eastern slant boundary and j = ^ for the western 
slant boundary. Table 2 shows the local coordinates in different regions. 

TABLE 2 
Local Coordinates in Different Regions 
Region Coordinates 






X 




1 


X 




2 


X 




5 


ll 


= £2 


k 


I2 


= £2 



z 



V^i 



^/^ [x 4 - .7(1 - z] , 



U = 61 Ml - z) 
-1/2 



Til = Z 
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It should be remarked that solutions of the flow field in region 1 and 
region 2 are not valid at the corners where the upper and lower Ekman layers 
merge with the slant side boundaries • As expected in the boundary layer prob- 
lem^ stretching the normal coordinates of their respective boundaries will 
cause difficulties at corners where two boundary effects are contradictory to 
each other « We notice that the flow in the Ekman extensions at the corner does 
not satisfy the slant side boundary conditions* There are four corners of 
dimension 0(ei/^x e|/^) where none of the above scaling is relevant. However^ as 
Greenspan and Howard (I965), Robinson (1959) and Hunter (1967) have found these 
corner regions to be small and relatively unimportant; the remainder of the 
flow can be solved without solutions of these corner regions. 

Analytical Approximate Solutions 

THE ZEROTH R-ORDER SOLUTIONS 

Temperature Distribution 

The zeroth R-order temperature can be obtained, from equation (20); which 
is simply the Laplace equation* The solution of equation (20) satisfying 
boundary conditions (2?)^ (28) is 

rr,(o) 1 r sinh ag ^ | /_^>, 

T^ ^ = 21^ " H ^^^ 

where a and H are abreviated constants specified in Table 5» The zeroth R- 
order temperature is purely conduction* Therefore it has no boundary layer 
effects. 



k^ 



TABLE 5 

Abbreviations for Constaxit Functions 

Abbreviated 
Notation 

a 

H 

G 
H2 
G2 

H3 
G3 
H4 

G4 
Ci 



Meaning 


Represented 


2n sTj 


sinh 2n 4y 


cosh 2n s/y 


sinh 2a 


cosh 2a 


sinh Ja 


cosh 5a 


sinh ij-a 


cosh Ua 


2Aa 


277 


-Aa 


277 


i|-TtAa 


27n/7 


-Sjt^Aa 


27H 


"kjiAo 


81n/^H 


jtAa 



27n/7H 



Velocity Field 

Interior Solutions 

In the interior region^ boundary effects are small. Equations (18) and 
(19)^ with the first terms neglected, become: 



46 



and 

° - = . 



dz 



which yield 



2 

and 






(o) 1 cosh az 

Y = - — sm 2tcx 

o ^ r H 



+ A^(x) , (58a) 



y^"^ = B (x) , (59a) 

o o 

where A (x)^ B (x) are integration constants and possibly functions of x* 
o o 

A (x) and B (x) can only be determined by matching up the interior solutions 
00 

to the solutions obtained from the bound^-ry regions* 

Boundary - Layer Solutions 

Region It— the upper boundary (free surface). In the upper boundary re- 
gion^ the stream function and the meridional velocity are: 

o ^ 

(o) (q) , (o) 

V = qV + iV 

and the vertical coordinate is stretched as shown in Table 2. Equation (52)^ 
using the local coordinates^ and to the order of £^ yields^ 

n6, (o) :n2,(o) 



^7 



Since, ^ir vanished exponentially away from the upper boundary, the "con- 
stants" of integration vanish. Therefore equation {kO) is simplified to equa- 
tion (3^). That is, in region 1, 



|f- + xt^°^ = . (Ul) 



The solutions of equation (kl) are: 






(te) 



where the m. are the roots of 



ra* + 1 = . (U5) 

These are 



["!]■ 



mi^2;»3^4 = ± exp ± i Ij- , {kk) 

Due to boundedness condition as ^i -> oo^ and since ^i, is always positive in 
the model, only the two roots with negative real parts are relevant: 

mi, 2 = - 1 ± i . (i^5) 

The other two appropriate boundary conditions are 



(o) ^ afi 

or equivalently. 



2.(o) 



^ = -v 2 — = at z = 1 



.2.(0) 



^ ° = ^ = at tx = . (46) 

All the C (x) in equation {k2) can be obtained, leaving only one constant 



Inherited frpm the interior undetermined. 

Equation (I9), to the order of e on the "boundary, leads to 



Since, iv is of significance only in the boundary region, the integra- 

(o) 
tion constant can be chosen as zero. After substituting i\|/ in and together 

■with boundary condition 



^ =.0 at ? = (U8) 

will yield: 

,v^°) = -eV^-^ sin 2«x(cos -^ - sin ■^) , (U9) 

^/2 v/i n/2 

and 

^^^°) = -ei«e"^^^''^sin 2«x cos ^ . (50) 

Region 2 — the lower boundary (rigid bottom)* After stretching the verti- 
cal coordinate, equation (3^4-) and another equation similar to (^7), satisfying 
all the lower boundary conditions 

,(o) . /o) , ^ . at C, - , 



we obtain 



^v(°^ = ey^>r2 Tte'^^^sinanx cos-^ , (52) 

n/2 



h9 



and 



^V = - ei ite*"^^/ ^sin 2jtx (cos ^ + sin ^) . (53) 

The integration constant in equation (52) has been absorbed into the integra- 
tion constant of the interior velocity^ ^n^"^^ ^^ satisfy the boundary condi- 
tion at z = 0. 

Matching the boundary solutions with that of the interior^ the boundary 
conditions in equations {h6) and (51) have also yielded the two integration 
constants in equations (27) and (28) as 

A (x) = ( ' — + e/ N/2Tt) sin2jtx , 

o r 

2 4yn 

and 

B (x) = ei Jt sin 2jtx- . 

Then the interior solutions in equations (38a) and (590-) becom^^ 



V 



(o) / COsh az - 1 1/2 7— X . ^ .ON 

= ( - £1^ \/2it) sm 2atx ^ (38) 



^ 2n/7h 

and 

t = £1 It sin 2jtx . (39) 

Region 3 and region 4 — the eastern and the western slant side boundaries. 
Using local normal component stretched coordinates in the side boundaries^ 
equation (32) to the order of e^ leads to 






50 



where j =^ 5 denotes the pastern boujadary^ J == ^^ the western boundary. Notice 

1/2 

that the boundary layer thickness is the order of 82 ^ where e^ = eK* K is 

constant as shown in equation (56) . 

Equation (19)^ to the order of e, also leads to 



,(0) 



S|v^°^ + >^y^ -1/2 M. 



(o) 



se! 



1 + \^y ^ 



^k 



= , 



(55) 



in region 3« 

Equations (54) and (55) satisfying the boundary conditions on the eastern 



slant boimdary, 



(o) _ J°) _ M. 



(o) 



= ^ 



U. 



= at I =0 



(56) 



and the decay conditioi:;^, noticing | o^s always negative in region 5, we have 
the solutions, 

"j..^ii|_^(l.e^x/^ cosi-) 



,(°) - 



1 - coshaTi, ^gV^^^ 



ll 



. sin 2JtX7(l - Til) , 



(57) 



and 



f = Bin sin 2n\y{l - \)e ^' (sin — ^ - cos -^ ) 



(58) 



Similarly^ in region k^ the solutions are 
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.(°) - 



cosh aiig " 1 
2\/7"h 



- e 



■/^ 



2« 



n/1c 



^ (1 + \^7) 



(1 - 



.e./^/^ 



COS — ^ ) 



• sin 2n\7(l - il„) ^ 



(59) 



and 



,(o) 



EiJt sin 2^X7(1 - 11 2)6 2 



^ ' (cos "-" 



(cos -^ + sin -^ ) 



(60) 



THE FIRST R-ORDER SOLUTIONS 

Temperature Distribution 

The first R-order temperature^ T ^ is governed by equation (23)^ with 
boundary conditions (29) and (50). The nonhomogeneous terms of equation (25) 
are the products of temperature fluxes and velocity components of the zeroth 
R-order solutions* Zeroth R-order velocities have boundary layer character- 
istics and these velocities are forcing the first R-order temperature « Though 
the zeroth R-order temperature has no boundary layer character^ it is expected, 
that at least part of the first R-order temperature will have boundary layer 
characteristics Besides^ the small viscous parameter £^ the symbol of bound- 
ary layer character^ is also coupled with the differentiated, highest order 
terms of the temperature in equation (25) • Therefore^ as indicated in equa- 
tion (51)^ let 



(1) (1) (1) 
T^ ^ - ^t' ^ + T^ ^ , J = 1.2,5A 
J 

where T is the overall part and .T is the boundary part of the first R- 



52 
order temperature ♦ 

The governing equation (25), then, can be "written in local coordinates in 
different regions as following: 

in the interior region 

_o (x) Qjt^ ^ ,a coshaz ^ ^x fr-,\ 

y2rp\ / ^ ..., — ^^^ 23TX ( — cos 2iTx -" 1); (61) 



in region 1, 



M) 



^i^ = - sx e ^' sin'=^ 2jtx smh a(l - eibi) 



S^ 



'i n/2H 



(62) 



(cos -^ + sin ^) + 0(e) , 



in region 2, 



&^ = ^l/^Kojt ,-W^/^ 3i„2 2„^ 3i„h (aej/^^^) sin ^ + 0(e) , (65) 
0^2 H ^ 



in region 5^ 



► / ; — 

aiT^^^ 1/2 aJt e ^ ^ , ^ . Ijl fn a cosh a% 
|ai ^ _ gi/2 _£L2 sin 2n\7(l - \) sin -i 1 - 



cos 1-!i{zJ^ ^1 + I - ^7 + ^7 Hi) + ^Y^ sinh aT]^ sin 2n (61+) 



in region h, 
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bfl^^ 1/2 ffit e ^/ . o ^ /-. -n N • I£ 



a cosh a.% 
H 



cos 211(62'^ Ig _ i + \y - xyT] ) + S^ sinh an sin 2rt (65) 
2 2 H 2 



1/2 1 

(£2 ^2 - - + ^7 - ^^riig) - 1 



+ 0(e) , 



The boundary conditions for the first R-order temperature in equation (29) and 
(30) require T^^^ = at z =^ 0, z = 1, and dT^^V^K = at x = ± 1/2. As we 
notice from equation (6I) and equations (62), {63), {6k), (65) that ^T is 
order of unity and .T^^'^'s, ^.re orders of e^'^. For simplicity, we can require 

J 

tha.t both T and T ^^ j = 1^2, vanish at the upper and lower boundaries 
to satisfy the boundary conditions -without losing good approximation. On the 
eastern and western 3ide^ since the boundaries are inclined and all the dom- 
inant forcing terms in equations (64) and (65) vanish on the boundaries ^ we 
require .T^^^ = 0, j =^ 3^ 4 at | = 0, instead of having ST^^V^I = in con- 
sistency with the basic temperature requirement on the eastern and western 
boundaries* 

Interior Solution 

The overall part of the first R-order temperature T can be obtained 
by solving equation (6l) subjected to boundary conditions^ 



T^^^ = at z = 0^1 ^ 
o 



^9^^'^ =0 at X = + i 
ax at X _ 2 



(66) 



5^ 
The particular solution of equation (6l) is 



T = T— cos 2itx + cosh az - — ^ ^^gj^ a.z cos knx * (67) 

^ ^ ^ i^ n/T'h 12 n/7h 



It is convenient to transfer particular solutions into new boundary conditions* 
For the homogeneous solution^ equations (6l) and (66) lead to: 

V^T^^^ = (68) 



h 



and 



rr.(l) m(l) 

T^ ^ = - T at z == ^ 
oh op 



T^'^ = - T^^^ at z = 1 , (69) 

Oh op 



;,t(-' 












h 










1. 




= 


at 


X 


= 


± - 


ax 










'd 



The solution of equation (68) satisfying boundary conditions (69) is 



T = - ' ' — + (1 - G)z - 7— cosh az cos 2jtx 



+ -^ cosh 2 az cos knx + . "* — ^ sinh az cos 2Ttx (70) 



^ aTc(G •- Gg) . ^ ^ , 

+ — A-, — fii^ smh 2 a:? cos 4itx j, 

12>/7'h H2 



where H^ G^, H^^ G^ are abbreviated constants listed in Table J* Therefore^ 
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the overall part of the first R-order temperature is 






— - [cosh az - 1 + (l « G)z] + -— 

47 



[ 



1 - cosh az + ^ * ' sipb az | . cos 2nx + 



■] ■ 



an 



izJy'n 



[■ 



(G - G ) 
cosh 2 az - cosh az + -^ ^ sinh 2 az 

^2 



cos hnx 



(71) 



Boundary Solutions 

The boundary layer equation;? of the first R-order temperature in differ- 
ent regions are the equations (62); (63)^ {6k), (65). The boundary conditions 

are .T »s vanishing in the interior^ and 
J 

.T =" ; d = 1^2^3^l4- on the respective boundaries. (T2) 

In region 1^ equation (62)^ after integrations and having satisfied bound- 
ary conditions^ yields 



= £1' sm 2:n:x 



1 - e 



-^ (cos — - - sin -^ 



• (75) 



In region 2, similarly^ ve obtain 



T^^) = 



(7^) 



Since the forcing terms all vanished at the lower boundary as shovm in equa- 
tion (62), equation (74) means T has no boundary layer modification near 
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the bottom as expected. 

In regions 3 and k, equations (6U) and (65) With boimdary conditions lead 
to 



^(1) ^ ^1/^ g^ 



"/2'(1 + A. ^7)7 



sin 2«X.7(1 - \) 



a cosh a^x 

H 



(75) 



cos 



2^^.7(1 - \) + ^^ sinh aili sin 211X7(1 - t1i) (e ^/ cos ":£ - l) , 
^ J n/2 



and 



^(1) 1/2 ore 



— " sin 2Tt\7(l - 'Hs) 1 + ~ cosh a^i^ cos 27(Kj 



[■ 



(76) 



(1 - ■^a) - ^17^ sinh ang sin 2«\7(l 

H 



^.)l(e-^-/^cos-^-l) . 



Velocity Field 

The first R-order velocity components are coupled with the first R-order 
temperature and the zeroth R-order velocities. From the governing equations 
(21)^ (22)^ (23)^ we obtain the equation for the first R-order boundary layer 



scaling as 



a^^, (O , ^ 



,2,(1) 



az' 



3k Sz ^ bx Sz 

.,(q) . (o) .,(0) (o) 
bz bx bx. bz 



bz 



b^^^°^ 



bz ^bz^ 



Sx d2= ' ^ ^ax^ az " Sx ax^dz ^ 



(77) 
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(o) (o) 

Noticing that the T has no Ipoundary layer character^ \|/ is order of 

(o) 
£ with boundary layer part having the same order^ ancj v is order of unity 

with boundary layer part of the order of £ ^ , the largest term on the right- 
hand side of equation (77) is estimated to be the order of unity in all the 
boundary regions. To th,e prder of £ in the boundary regions, equation (77) re- 
duces to the same boundary layer scaling equation for the zeroth R-order, that 
is, equation (JS)* Thus, we can solve the velocity field of the first R-order 
in exactly the same way as "we did for the zeroth R-order solutions. 

Interior Region 

In the interior region, since b \|/ /Sz = 0, all the zeroth R-order terms 
in equation (2l) are zeros. Therefore, in the interior region, equation (2l) 
becomes 

o __ o 
which together with T in equation (7l) yields 



(i) a fsinh az ^ (G - l) cosh az 
v' ' == 7~" — — — - 2:itz - -^ '^— ' — 



sm 2nx - 7-— 
D7H 



(78a) 



sinh 2az . ^ ^ (G - Gg) ^ ^ 

- smh az ^ • — — — ^^ cosh 2az 

2 2H2 



Equation (22), to the order of £, leads to 



sin hfix + Ai(x) 



:. (1) :. (o) :. (o) 

Q _ O^ O 

Sz bx bz 
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which, -with equatipns (58), (59)^ yields 

(1) Tt^ 
\|; = £1 "~ cosh B.Z sin knx + Bi(x) . (79a) 

The integration constants will be determined after matching the solutions. 

Boundary Regions 

In region 1 using the same local coordinates as shown in Table 2, equa- 
tion (77), to the order of 8, can be reduced in similar manner to the similar 
equation (5^) in the zeroth R-order solutions « With the boundary conditions 
similar to (46) and (48), equations (77) and (22), to the order of e, lead to 

^v^^^ - , (80) 

and 

^^^^^ = . (81) 

The velocity component due to ^T modification to the overall temperature 
distribution on the bounda-ry has been neglected within the approximation ac- 
curacy. 

In region 2, similarly, w^ obtain. 



(^) 1/2 jt^d - G) -ia/ ^ . , ia /Qo^ 

V = ei — ^^ — ^ e ^^ sm 4tcx cos -^ (82; 

and 

(1) Tt^(l - G) -Wn/2" . , , ia ^ . i^^ fQ^\ 

2^ = - £1 — "^'~' ^ e "^^ sm knx (cos -^ + sm --fj . (83) 

2 n/t^H n/T n/2 
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After matching the interior ^.nfl the boimdary solutions on the boundaries^ the 
two integration constants, A;i^(x) and Bi(x) have been determined. Equation 
(78a) becomes 



.(i) 



a Fsinh az ^ . (G - l) ., , v 

= r- ^ 2itz + -^ -^ (1 - cosh az) 



p 



^^ ^/7 



n/7 H 



sin 2Ttx 



ML 



sin h 2az . ^ ^ (G - G^) . , ^ .1 
- si,nh az + — '"^' (cosh 2az - 1) 

d 2^2 J 



(78) 



g^ . .. >M — ^ i^ Sin 4JCX 
\i2y H 



Equation (79a) becomes 



(1) jrfsinijjt^ . ^ ^\ 

f = £1 ^ (cosh az - G) 

2n/7"h 



(79) 



In regions 5 and h^ noting that the order of unity part of the zeroth R- 
order meridional boundary velocity, v and v , are functions of the un- 
stretched tangential coordinate rj only, equation (77) leads to an equation 
similar to (32) as mentioned before* Then, in similar manner, we obtain, in 
region 5 



^-^ = , gV^^ ^ e ___ ^^^ 4Tr\7(l-Tli) cos — (cosh aTli - G) -h{\) {8k) 

\f2jE{l + X^y) sf2 



where 



6o 



M.) - ty 



sinh a**! ^ (G - l) , J] 

' - aixn + -^ ^ (1 - cosh a-q) | sin 2it\7(1 - y\) 

V 7 vy H 






L - ^ -^ n/27 H 



62 — — (cosh ari - G) }>sin ^^11X7(1 - tj) ; (85) 

n/27'h(1 + X^y) 



and 



(1) n^e ' 



= ei 



sin h^yiX - \) (cosh aTJi - g) (cos -^ - sin -^) ,(86) 



2\/7"h n/2 ^2 

in region 4, 

^v^^"^ = £2^^ -^^ ^ ^ sin iin\7(l - rj^) cos -^ (cosh a^l^ - G) + hi'n^) 

n/27h(1 + X^j) v/F ^Q^^ 

where h('n2) is a direct substitution of "02 for r\ in equation (85) and 

/ \ 2 "2/ V 2 ? & 

A|; = - £1 -^^^ — ^ ' sin knXjil - %) (cosh aTj^ - G) (cos -^ + sin -^ )^ 

2n/7"h v/i" \/F 

(88) 

1/2 
Note that all boundary layev solutions in our model study are weaker by 0(£ ^ ) 

than those obtained in the case of rigid surface (see, for example, Robinson, 

1959) • This agrees with Hide's (196^1-) result that the Ekman boundary layer 

associated with a free surface is weaker than that with rigid surface. 
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Estimated Errors of the Approximate Solutions 

TEMPERATURE 

(2) 
The order of magnitude of the second R-prder temperature^ T ^ can be 

evaluated from its governing equation (26). Frpm the zeroth R-order solu- 
tions^ we know i is 0(e) both in the interior and boundary regions; v 
is 0(1) in the interior and 0(e ' ) for the boundary layer part; and T 
has no boundary layer characteristics* The first R-order solutions have ex- 
actly the same order of magnitude as that of the zeroth R-order ♦ Therefore^, 

i/a 
equation (26)^ having the forcing terms 0(e) in the interior and 0(e ' ) on 

the boundaries^ will lead to the same results as the first R-qrder temperature^ 

(2) (2) 1/2 

namely^ T^ ^ is Q(l) and .T^ ^'s 9.re 0(e ' )• Hence the highest order of 
o J 

(2) 
T is 0(1). Thus^ the validity of th^ series expansion of the temperature 

in powers of R as shown in equation (17) has been verified. The asymptotic 

character of the infinity seriejs is now plausible. 

To the expansion of second power^ the approximate error of temperature is 

O(a^R^) which is rather small, 

VELOCITY 

The meridional velocity is dominant in the flow field- It is 0(l) in 
both the zeroth R-order and the first R-order* Therefore^ the approximate 
error for v^ to the second power of the infinite series^ is 0(R^)» A similar 
approximation error is found for the stream function. 
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Type -A Circulation 

Type-A circulation is expected to occur during the summer -heating period 
when the temperature in the lake has already risen above 4 °C with the middle 
part of the lake still colder than the edges, and during the winter-cooling 
period when the lake has cooled down to less than k^'C and is still cooling. 
The solutions of the flow field have been obtained by applying the surface 
temperature of the lake a$ a concave cosine function as mentioned in the last 
chapter. Notice that the thermally induced force during the winter -cooling 
period, since the fluid is less than k°C and the air temperature is cooling 
the lake, has exactly the same effects as the summer -heating period, except 
the former has smaller AT and in the negative sense. 

Type -B Circulation 

During the spring -heating period and the autumn-cooling period the sur- 
face temperatures of the lake and their induced thermal body forces have ex- 
actly the reverse forms and opposite effects as that of the winter- cooling and 
summer -heating periods in the type-A circulation. Therefore, type-B circula- 
tion is the reverse form of type-A circulation. 



IV. THE THERMAL BAR CURRENT 



Mathematical Formulation 



The general governing equations (l), (2)^ (5) with boundary conditions {k) 
and (6) are the differential system for the flow field during thermal bar periods 
of the Lake Michigan model. The current circulation pattern induced by the ther- 
mal bar phenomenon specified in Chapter I is called the thermal bar current or 
the type-C circulation in our terms. The thermal bar phenomenon is common to 
all large dimictic lakes during the time of spring warming-up period or the 
autumn cooling-off period. The latter^ due to the smaller horizontal tempera- 
ture difference, is much weaker. Since the effects of the thermal body forces 
in inducing currents during these two pe^riods are exactly the same, we shall 
consider the spring thermal bar phenomenon as representative, 

THE APPROXIMATE EQUATION OF STATE 

During the prevailing period of the thermal bar phenomenon, the tempera- 
ture near the shore of the lake is above k°C while the temperature in the mid- 
dle portion of the lake is still less than the temperature of maximum density* 
The relationship between pure water density and temperature is shown in Figure 
8. Generally, Lake Michigan is rather shallow. The deviation of the maximimi 
density of fresh water away from 4°C due to the hydrostatic pressure can be 
neglected. The density of the natural fresh water in the lake follows the 
curve closely as temperature changes (Welch, 1952). Therefore, within the ther- 
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Figure 8. Relationship between freshwater density and temperature (after 
Welch, 1952). 

ml bar temperature range (from 0°C to 12°C, for instance), the density vari- 
ation with temperature change can be approximated as a cosine function. Let 
us define 



p = p cos A(T' - T ) 
^ o o 



(89) 



where p is the density of water at any .local temperature, T', p^ is the density 
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of water at reference temperature, T = U°C, and A is a constant which adjusts 
the equation of state (89) to the nearest good approximation of the density- 
variations within the temperature range (0°C < T' < 12°C). The appropriate 
constant A is^ in general, a small n-umber (10""'^ to 10"^). Thus the rate of 
change of density with respect to temperature is a function negatively corre- 
lated with the absolute value of the temperature difference between the local 
temperature and the temperature of maximum density. Then the thermal body 
force is 

— = cos A(T' - T) - 1 , (90) 

^o 

which satisfies the special requirement that the water density decreases as 

temperature varies away from ^°C. This approximation is used in the type-C 

circulation. 

NOKDIMENSIONAL GOVERNING EQUATIONS 

These follow closely the processes in Chapter III, using the same non- 
dimensional variables as equation (9) and with the same basic assinnptions as 
equations (U) and (12), the nondimensional governing equations of the thermal 
bar current, after introducing the nondimensional stream function^ become, 

(91) 
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where ^, V^ are operators and jy e, a are nondimensional parameters (Table 2). 
Note that the basic geostrophic thermal wind relationship is now a nonlinear 
equation, 

2fl -^ = Ag sin A(T' - T^) "^ , (9^^) 

which after being nondimensionalized becomes 

^^ V Sv AgAT sin AAT . _ ST 

where the sin AT has been normalized in order to make the two terms of the same 
order in the nondimensional equation. Thus the dimensional reference velocity 
becomes 

and the thermal Rossby number is now 

V_ _ AgDAT sin AAT . ^. 

Similarly, since R « 1, by observation, the regular perturbation method can be 
used to achieve approximate solutions. All dependent variables are expanded 
in powers of R as shown in equation ( I7) . The same argument in Chapter III can 
be used here for the validity of expressing the temperature in powers of R, We 
also expand sin AT in a power series. Due to the smallness of A, terms higher 
than second order in the sine series have been neglected. Then, after equating 
each power of R in equations (91), (92), and (95)^ the zeroth R-order equations 
are: 
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^(o) ^^°^ .Jo) ar^°^ ■ ,^^. 

,va^(o) _|/!l = o , (98) 

dz 

v^r^^ =0 . (99) 



We notice that from equation (99)^ the zeroth R-order temperature is still due 
to conduction only. 

The first R~ order equations are: 

>. ( O) -n3 ( O) -N ( O) n3 ( o) 

dz c3z dx dx dz 

and J 

/ \ ^ ( o) ^ f o) ^ ( o) ^^( o) 

dz dx dx dz 



The second R- order equations are: 

'^'^ 3z ^ ^ Sx Sx Sz SxSz 



a^.^iill ^l!!^iLl ai^^>k^°^ .A-^i^ 



(o) ^3 (i) 
Sx dz^ ' Sz SxSz^ 5x Sz^ ' ^ ^ Sz ^^ 






^^ir) . (103) 
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(lOlt) 



^Sz ^ 3k Sz Sz ^ Sk Sz 



(105) 



Boundary conditions, with reference to equations (k) and (6), are: 

T^°^ =i(l- 2cos2rtx) . t^°^ =f^^=^ = at z=l , 

T^°^ =0 , and Mr^°^ = v^°^ = ^ = at z = , (106) 

and q = on the slant boundaries and 

for the zeroth R- order equations; 
and 

(i) (i) £i^ ^^^^ , , 

T^^^=0 , t^^^ =|^ = v(^^ =0 at z=0 (107) 

^ "^ dz 



.(i) ^-^^ 1 

q = on the slant boundaries and -r- = at x = ± — 



->i 



for the ith R-order equations, i = 1,2,5,... 

Equations (97), (lOO), and (IO5) are highly nonlinear in temperature, 
which produces the thermal body force in generating currents. Fortunately, by 
virtue of the energy equation, the temperature field of every ith R-order equa- 



69 

tion can be solved first and the nonlinearlty in the ith R- order momentum equa- 
tion is eliminated. In this way^ the temperature and velocity fields can be 
found. 

BOUNDARY LAYER SCALING 

Equations (97)/ (98) with the Ekman number £ « 1^ lead automatically to 
the method of matched asymptotic expansion as before. Using the same notation 
to indicate the ith R-order dependent variables in the interior as well as in 
the boundary regions (Figure 7) as shown in equation (51)^ equations (97)^ (98)> 
and (99) result 

(o) 

Since T has no boundary layer characteristics as we noticed from equation 

(99) which is simply the Laplace equation, the right-hand side of equation (108) 
is 0(£) throughout the flow field, even in the boundary layer regions. After 
stretching the normal coordinate in each boundary region, the right-hand side 
of equation (lOB) is relatively only 0(£^) which can be easily neglected without 
causing much error. Then the characteristic equation for boundary layer scaling 
is exactly the same equation as (52). Following the same steps, we obtain 
the boundary layer thickness in region 1 and region 2 is 0(ei/^) and in region 
5 and region k is 0(62/^). All local coordinates are as shown in Table 2, 
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Analytical Approximate Solutions 

THE ZEROTH REORDER SOLUTIONS 

Temperature Distribution 

The zeroth R-order temperature solution satisfying the Laplace equation 
and the boundary conditions (106) Is: 

T^ = -(z ^ - sinh az cos 2Tfx) (109) 

which has no boundary layer characteristics. 

Velocity Field 

Following closely the procedures of the last chapter^ interior solutions 
can be obtained by neglecting the viscous effects. Equations (97) and (98) 
in the interior region^ become 

and 

■i — = ° • <^^i' 

Since T is known^ both equations (llO) and (ill) can be integrated. The 
integration constants will be determined after matching the interior solutions 
with the boundary layer solutions. As before^ we are using the technique of 
additive type botindary layer analysis. In the boundary region^ after stretch- 
ing the normal coordinate^ the governing equations in terms of the boundary 
local coordinates can be treated as ordinary differential equations. In each 
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boundary region^ the stream function can be obtained first from the equation 
similar to equation {^k) by satisfying the decay condition and the boundary 
conditions* Then the meridional velocity in the boundary region can be found 
from equation (98) after substitution of the boundary layer stream function and 
satisfaction of the meridional boundary condition on the boundary. Considering 
that the interior solution is dominant in the interior region and has no bound- 
ary layer characteristics and that the boundary layer solution is important only 
at its respective boundary^ the zeroth R-order solutions have been obtained as: 
in the interior 

( o) ^TcA z J 

V = — t-CCtT" cosh az - — ^ sinh az) sin 2jrx 
o 9 ''Ha Ear ^ 



(TS^T sinh 2az - -A-) sin kjcx] - e^/^ ^f2 b (x) , (112) 



^H^a ""^^^^ """ " 2H^ 



where 



in region 1^ 



y°^ = eibjx) , (115) 



^TtA 

b (x) = — — ( sin 2Trx - sin kitx) ; (H^) 



,,(°).,iA!e^e-^./Vi(,,„ii.,,,ii, ^ (,„) 



in region 2^ 



V2 n/2 ^2 



it^°^ =-eib (x) e-^^/^cos-^ ; (ll6) 

° 42 



/°) = s.V-^/ib (x) e-^-/^ cos -§.. , (117) 

o st2 



72 



,^(°) = - eib (x) e'^^/"^ (cos -^ + sin -^) ; (ll8) 

o V2 n/2 



in region 3^ 



3v(°) = .^M . eV= ^ .^,,,,[^ .^|^(ee./-/i ,„, i. . ,,j , 






(119) 



where 



3t^°^ = eib (ni)|^/^ (sin ^ - cos ^) ; (120) 

° V2 ^/2 



^ (n) =-7r[sin i^jr\7(l - n) + sin 27t\7(l - t]) ] ^ (121) 

o 9 



and 



_ , V UjtA ^^ sinh ari ti cosh ari x . o ^ / •, \ 

in region k^ 

J^) =-d(^) -eV-^/ib(Tl.)[i+■r-^^(e-^-/^cos-^-l)] , 

(123) 

4t^°^ = eib {^2)e^^^ (cos -^ + sin ■^) (12l+) 

° ^ v/2 

where b (t]^) and d (rj^) are direct substitutions of t]2 f*or t) in equations (121)^ 
o o 

(122). 
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THE FIRST R-ORDER SOLUTIONS 

Temperature Distribution 

The governing equation and boundary conditions for the first R-order tem- 
perature are the same as for the general thermal currents. We notice from 
equation (102) that the first R-order te^mperature has boimdary layer character- 
istics because of the forcing terms due to heat fluxes coupled with the boundary 
layer velocity components. As before^ sifter stretching the normal coordinate 
and rewriting in terms of the local coordinates as shown in Table 2, the bound- 
ary layer part of the temperature can be obtained by integrating equation (102) 
with all the boundary layer effects included and by satisfying boundary condi- 
tions in (107) on its own boundary. Using the same notions as indicated in 

equation (jl)^ the .T *s are summarized as: 

J 

in region 1, 



^(1) 1/2 2aaT . / X . 

T' ^ = £-,/ -r-- b (x) sin 2Ttx 
1 -^ ^ o ^ 



'^1/^ 



1 - e 



(cos -7^ - sin -^ ) 



n/2 



^/F 



> (125) 



in region 2, 



(1) 



(126) 



in region 3^ 



knkj 1 ^i/^ I 

+ —^ sinh ar)^ sin 211X7(1 - ti^) (e cos — '" - l) , 



(127) 



7h 

in region k, 

(i) I ^'^V^^ 2a 

(128) 
- -- sinh a'n2 sin 2%\j{\ - *n2)](e ^^^ cos -^ - l) , 

where b (x) is shown in equation (llU)^ b (tj)^ in equation (l2l). In the in- 
terior region^ after substitutions of equations (IO9) and (ll^i-)^ equation (102) 
leads to: 

^ (1) STT^AcT/a . 2a 

V T = — — — ("-; cosh az - cos 2jtx + 2 cos 4Ttx - — ■ cosh az cos 2jtx 

(129) 

+ — cosh az cos ^irx - — '- cosh az cos 67rx) 
H H 

The boundary conditions are still as expressed in equation (I07). Since the 
homogeneous solution of equation (129) has terms like 

cosh naz cos 2niTX ^ n - 1^2^5,...^ 

the particular solution^ knowing that equation (129) is an even function, is 

T = Ci cos 2Ttx + Cp cos ^TTx + C3 cosh az + C4 sinh az cos 2tcx 
o p 

+ C5 cosh az cos ^ttx + Ce cosh az cos Sttx . (I50) 

Substituting equation (I50) into (I29), we obtain 



2Aa 
''^ " 277 ' 


^^ " 277 ' 


27n/7 


-8^Aa 
C4 - J 


-i|-rtACT 


«A0 



(151) 



75 



For convenience^ all the C. *s are listed in Table 3 as abbreviated con- 
stants for later use. 

After the transfer of the particizlar solution as new boundary conditions 
for the homogeneous solution^ we obtain: 

T = 03(1 - G)z - C3 - Ci cosh az cos 2Ttx - (C^ + C5) cosh 2az cos hitx 

c* c* " n 

- Cg cosh 5az cos 6:rtx + {— - C4) sinh az cos 2Trx 

H 

(Gg + Cc,)G ■ Cg - Cc,G . ^ ^ , ^ C6(l - G) . ^ ^ . 

+ -^-^ ^^ ^ ^^ smh 2az cos ^ttx + ^^ ^ smh 5az cos bitx 

H2 H3 

(152) 

Therefore, the overall part of the first R- order temperature without boundary- 
layer characteristics is, 

T^^ = 03(1 - G)z + C3(cosh az - 1) + {Ci(l - cosh az) - [04(1 - z) 



„ _iA Z.-| sinh az] cos 2Ttx + [C^ + C5 cosh az - (C^ + C5) cosh 2az 

(C + C ^C - G - G G 
+ -^— ^ ^^ ^ ^ ^— sinh 2az"] cos 4irx + [Ge(cosh az - cosh 5az) 

+ ^^ ^ " — ^ sinh 5az] cos 6itx (I55) 

H3 

and has 0{e^'^) modifications on the boundary regions. 



Velocity Field 

In the interior region, since S ^ /Sz = 0, equations (lOO) and (lOl) lead 
to 
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¥^ = ^^'^'#^^'°'#. > 



and 



:. (i) ^ (o) . (o) 

o oj^ o 

Sz Sx Sz 



(135) 



The values of T ^ v y V 9 and r are known as shown in equations 

000 

(27), (58), (59)^ and (155)- Therefore the solutions for v^^^ and J"^^ can 

o o 

be obtained by integration. As before^ the integration constants will be de- 
termined after matching the boundary solutions and the interior solutions. 

In the boundary regions, equations (lOO), (lOl), and (102) lead to the 
boundary layer characteristic equation: 

Sx Sz Sx Sz ^ Sx Sz ^ 

^ ^^ Sk ^^ 7 Sz SzSk 
~w,(o) n2,^(i) -n (o) -^3 (o) -^ (o) .3 (o) 

Sz azSx ^^ ^"^^az axSz^ " sk sz^ 






^ ( aT- a^^ - ■^- ^fe-) ^ • (156) 

The right-hand side of equation (I56) is separately denoted as I, II, III, IV, 
and V for each successive part for convenience. It is easy to check from so- 
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lutions obtained a*bove that I is 0{e^/^), II is 0{e^/^), III is 0{e^/^), IV is 
0(e), and V is OCeV^) in regions 1 and 2. In i^egions 5 and U, note that the 
contributions normal to the slant "boundary of each individual term in the dif- 
ferent parts, except part III, cancel each other and what is left is a portion, 
0(e^/2) lower than the normally differentiated portion. It is also noticed 
that the portion of order of unity of the zeroth R-order "boundary meridional 
velocity, aJ \ is a function of r| alone. Thus every part on the right-hand 
side of the equation (I36) is estimated to "be 0(£V^). Hence, to the 0(e^/^) 
on the boundaries, equation (I36) leads to the same boundary layer thickness 
and same governing equations as before. 

Following closely the procedures in the zeroth R-order solutions, after 
having matched the boundary solutions to the interior solutions and satisfied 
all boundary conditions, the first R-order solutions for the velocity field 
are: 
in the interior region, 

(1) 2TtA,Ci 2C4 18C.^ + %Cp + 17Cc, „ P 2C^(1 - G) - Ci(G - l) + C^H 
o^ =— {^^--^ -H ^^^^ -C,z - [ ^^3 

, (Cp^C.)Gp-Cp.C,G ^.^^ ^^ ^ iC,^C.)^,^C,^C^^ ^.^^ ^^^ 

^a.^ a 2Ha ^ UHa 

(Cg + Cc,) . ^ ,-^ 2Ci, . , 2Cc,(l - G) -G^fG - l) +C4H 
- ^ "^^ — ^ cosh 3az + (— ^ + — ^) z sinh az + — ^ ^T^ 



6Ha ^ a a Ha 

-^ z cosh azj sm 2itx - — "( — ^ — ~ TTTs 
a 3 ^ ^Ha 



+ 2^C; - ^P + 2C^(G - 1) - 2C4H ^ ^ ,c^ ^ c^ ^ i^ w£ 
^ 8Ha H^ ^ "^ ^ H .U ^^^''^ 

_ Cp + C,G - (C^ + C,)Gg _ C^ ^ Ci(G - 1) - C4H _ MG3- G) ^.^^ ^az 
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- ^it^^ -^^ ^^- - (V^ +%) eosh az -. (-^ +^a4^ 

SHHga '•Ha a"^ '^ ^li-Ha a^ 

, ^(Cs ^ C.) ^.^^ ^ ^ 2C.(G, - 1) . 2C.(ap - G) _ C^ 

a ' ^ Haa 2Ha-' •' 

. 3in l.«x - i^{Q^ -^ 1^^, -^ TC. ^ 3(C^ ^ C.)Gp - 3Cp - ?C.G ^.^^ 
5H 59- ^a 2H2a 

2H2a JHsa -^ ^ b^ 2a ^ 

• cosh az - -7--^ cosh 2az + (-^^— — ^ + — -|) cosh 5az + -^Si^ z slnh az 
^a ^a 5a a 

- ^ z sinh 5az + "^^j^^^ - ^^ z cosh 5az3 sin 6nx - ^ M^ ' ^3) 3^^ gaz 

a nsa ^H oHsa 

+ ^^^ow I si"^ ^^2 + ^ cosh 2az - -^ cosh kaz - 7^] sin 8jtx 
12H3a 5a 12a 4a 

- el/^ 42 bi(x) , (157) 
where 

f . l6iT^A^ 2G 2 H2 ^ 1 . . , ^,G Ix . , ,2G 2 

G - 1 1 H 

"^ "iFT " "i?^ ^^^ ^"^^ ^' ^2?I " ^^ ^^^ ^"^^^ 5 ^^58) 

and 

f 1) l6Tr^A^ ,,2, a - z cosh az sinh az - H. sinh 2az - Ho 1 - z, 

o^ = ^- -sT-f^i^ 1 ^ ? ) ^ Wl — ^ ^F"^ 

. 3in 2:tx - ( ^ "°^\^^ - ^ - ^^"^^;f - ^ ) sin 4^ -h t|(^-COsh_a^^ 

sinh az - H , , cosh az - G z-l_ . . ,sinh 2az - Ho , 

- -^ ) + — ^ga; - -^^] sxn 6rtx - ( 21^^ - z + l) 

• sin 8x} ; (159) 
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in region 1, 



iv^""^ = , (140) 



11^^"^ =0 ; (lIH) 



in region 2^ 



2V^^^ = 42 eiA b^(x) e"^^/"^ (cos % + sin ■^) , (ll^2) 

n/? n/2 






2^l^ = -£ibi(x) e ^ cos 



in region 3^ 



-"^'^ = ^-'^' TITI^ b3(ri.)(e^-/^ cos i^ - 1) -H d3(T,x) , {±hk) 

si^^ = eibaCrii) e^^/"^ (cos -^ - sin % ; (li^5) 

N^ n/2 



where 



■u ( \ ^ -l6?t^A^ 2. sinh an - H ti cosh an - G v sinh 2an - Hg 1 - n -i 

o ^ /T \ , / Sinh an - H n cosh an - G» . , , ,, . 
. sxn 2rt\7(l - n) + ( :^ - -" :^ ) sm I^jrX7(l - n) 

, ^2^ sinh an - H n cosh an - G^ G - cosh an n - 1-. . /- /. \ 
+ [-( ^ - ^ ) + -^^^ + ^^] sxn 6«X7(1 - n) 

+ ( ''"gH^;^ " "^ - T, -H 1) sin 8«Xr(l - ri) , (l46) 

. f^^ =^2A,c,.2 + 203(1 - G) - Ci(G - 1) + C4H ^ (Cg + C^)Gp - Cp - C^G 

(Cg + Cc,)Gp - Cp - CsG . ^' ,Ci J+C3_+9C2+3C5 2C4, 

• smh an - ■^— "^ ^^ ^ sinh Jan + (2 + — ^ " ^ "^-^ + —^) 

' bHHpa -^ ' ''a^ 2Ha a^ '' 



20. 

2 
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. (cosh ari - 1) - -^^^^Ccosh 2aTi " l) - ^^^^^^^ (cosh Jar] - l) - (^^ 

^ 2Civ . , 2Cq(l - G) - Ci(G - 1) + GaB. Ca o 

+ 1~)ti sxnh ail - — ^^- ' ^J-^ i ^^ T) cosh ari + ^ tj^ cosh ar,} 

• o % /-T \ itA ,2Ci(G - 1) - 2CaH ,C4 Ca , x a 

. sm 2rt\7(l - T)) - —[ ^ ^ ^ ^ + (^ +-A + i^C2)Ti2 

■ Axm — ^ sinh 4aTi - (— — + -^)(cosh ai] - l) + {-rrf^ + "a — ^ + — \ ^ ) 
OHHga ' ' Ha a'=^ '' ' ^UHa'^ a'^ 2Ha ^ 

■ (cosh2aTi-l) + ^(cosh i^a^l - 1) + ^ t) sinh ari - [77^ + ^(^S + G^) 1 

oha a 2Ha a 

• sinh 2aTi + [^^g^^P - ^^jj^/^"^^^ ' ^^ - 1^] Ti cosh 2aTi3 sin l^^Xrd - ti) 

_ 2«A,3C^ + ?C.G > ?(Cp H- C^)Gp .(C^ + COGp - Cp - C.G 

3H ^ 2H2a ^^""^ ^"1 •- 2H2a 

^ ^^^^^^ ^^"^ 5aTi + (^ -f ^22^^)(cosh an - 1) + ^(cosh 2a, - l) 

/ Cg + Cc, HCew , , ,x 3HC6 , ^ ^ HCs . ^ , 

- {-^ + T2)(cosh Jar) - 1) - -^^ — ^ smh ari + — ^ sinh Jari 

^a ja a a 

HCsCGq - G) u , 1 • /■ X /-, N rtA^CsfG-Ga) . , ^ 

■ — ^fZ ^ ^°^^ ^a,} sxn 6jt\7(l - t]) - — { ^V„ ^^^ sinh 2aTi 

■tlsa en 0113 a 

. sin 8nX7(l - Ti) + sl/^ 8I ^^ " F " te^ '' SF^ ^i" ^nXyd - ,) 

/^ 1\ . 1 X ^n \ /2G 2 G - 1 1 , . , , , 

- ^i; - ^^ sxn WA.7(1 - n) - (- - ^ + ^p^ - ^) sxn 6«Ml - ,) 

- (^ - 1) sin 8rt\r(l - Ti)] ; (147) 
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and in region hf 

^^^'^ ^ ^^ (ifx^^) ^3(^P)(1 - e-^^^ cos ^ ^ deCria)) (1^1-8) 

^^(^) = ^SibsCris) e^^'^ (cos -^ + sin •^) (l49) 

all b3(r|i 2)^ ^aC'Hi 2) ^-^^ direct substitutions of t)i 2 f'o^ T i^ equations (1^5) 
and (1^5). 

The magnitudes of approximation errors for the higher order terms neglected 
in the power series expansion in terms of R can be estimated from the second 
Reorder equations (105)^ (104), and (I05). They are O(cr^R^) for the temperature 
and O(R^) for the velocity field. 



V. EDDY VISCOSITY AED EDDY DIFFUSIVITY IN lAKE MICHIGAN 

Introduction to the Characteristics of Eddy Viscosity 
and Eddy Diffusivity 

Theoretically, a flow field of incompressible Newtonian fluid can be found 
by solving the equation of conservation of mass and the Navier-Stokes equations 
of motion, and satisfying its respective boundary and/or initial conditions. 
In the event that a transformation of energy is important, the energy equation 
must also be used. When the flow is viscous and laminar, there appears in the 
equations of motion a term related to shear stress which, according to Stokes' 
law of viscosity, is 

du 
'^^^ ^^ 

where the subscript 1 denotes "laminar," z is normal to the direction of the 

flow velocity u, and [x is the molecular viscosity which is a property of the 

fluid depending only on temperature under ordinary conditions. In compressible 

laminar flow, when temperature variations are involved, the heat flux following 

Fourier's law of thermal conduction is 

dT 

Qi = -k — 
dz 

where k is the molecular thermal conductivity which is also a property of the 

fluid and usually depends only on temperature. 

On the other hand, if the flow is turbulent a kind of motion with random 

irregular fluctuations is superimposed on the main flow. Since its fluctua- 
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tions are irregular, it is impossibly to describe the motion in all details 
as a function of time and space. However, one can study its average behavior 
using the statistical methods which have been pointed out by Hinze (1959): 
''Turbulent fluid motion is an irregular condition of flow in which the various 
quantities show a random variation with time and space coordinates, so that 
statistically distinct average values can be discerned." Therefore, in de- 
scribing a turbulent flow in mathematical terms, it is convenient to separate 

it into a mean motion and a fluctuation, such as 

->• "^ ->• 
q = q + q' , 

(150) 

T =; T + T' 
where the "bar" indicates mean quantities and the "prime" indicates the fluc- 
tuation from the mean. 

The presence of random irregular fluctuations often manifest themselves 
as an apparent increase of viscosity and diffusivity (Defant, I961). In tur- 
bulent flow, therefore, additional Reynolds stresses must be considered in the 
general governing equations of the motion. 

Analogous to the coefficient of molecular viscosity, Boussinesq (I877) 
introduced a. mixing coef ficiei.t , A for the Reynolds stress in turbulence to 
express it in terms of the mean velocity gradients in the flow field 



— — — ^ du 
x^ = -p u w = A^ - 

where the subscript t denotes turbulent flow. In a similar definition to that 
for molecular diffusivity, Schlichting (1968) postulated that the turbulent 
heat flux is 



&k 



dT 

Q. = -C A fi 

t p q dz 



The exchange coefficients of turbulent momentum and turbulent heat flux^ A 
and A ^ both have the dimension of a viscosity, \i, (g/cm/sec). 

If the turbulence is isotropic (i.e., the relation between stresses and 
strains is independent of choice of axis) then A corresponds to the molecular 

"G 

k ^t 

viscosity jj., A corresponds to the molecular heat diffusivity — ; and v - — ^ 
^ Cp t p ^ 

the eddy kinematic viscosity^ corresponds to the molecular kinematic viscosity 

V = — . Then the turbulent shear stress is 

I P 

du 

T. = P V , — . 

t t dz 
The eddy kinematic viscosity, v , can be expected to occur as a constant only 
if the turbulence field is homogeneous. Generally speaking, since the eddy 
viscosity is 10^ to 10^ greater in magnitude than molecular viscosity, the 
latter may be neglected in many cases. As pointed out by Lamb (1932) in illus- 
trating the wind producing current, in reality jj. is replaced by A , 

Using the analogies mentioned above, it is possible to attack geophysical 
fluid problems which are turbulent in nature. In this way we can probably 
predict or evaluate the time-mean value of the entities of these complicated 
turbulences whose complete theoretical formulations have so far proved impos- 
sible (Lamb, 1932 ). 

Nevertheless, not all the turbulence is homogeneous and isotropic. The 
eddy viscosity usually does have directional preference. Eddy viscosity be- 
comes effective only if there is some flow in the fluid and its value depends 
on the velocity distributions and on the characteristic length of the flow 
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field as shown in Boussinesq's hypothesis and von Karman's similarity theory 
that 



A «2 l^^l -P |dU| 

K = p ^ \:r\ , or V, ^ r — 

t 'dz' ^ t 'dz 
where i is the Prandtl's mixing length. Therefore^ A or v. is not a property 
of the fluid itself, whereas the molecular viscosity, |u, is. However, it is 
still extremely useful to introduce the eddy viscosity in relating the Reynolds 
stress to the mean velocity gradients of the flow. As A or v is not a phys- 
ical constant, but varies from one part of the flow to another, estimating the 
magnitude of the average eddy viscosity is the key factor to the success of 
geophysical fluid modeling. In Lake Michigan, though the mean velocity is in 
general less than 1/2 m/s (Ayers, et al., 1958), the flows are turbulent in 
nature. Therefore, it is important to estimate the order of magnitude of the 
eddy viscosity and eddy diffusivity. 

Theoretical Background of the Formulae 

In order to evaluate the magnitude of the eddy viscosity and eddy diffu- 
sivity, different approaches have been used for the easier management of 
appropriate field data . 

REYNOLDS STRESS AND VELOCITY GRADIENT 

Analogous to the kinematic molecular viscosity formally introduced by 
Boussinesq (1877), Reynolds stresses are written, in the conventional Cartesian 
tensor notations. 
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T. . = -P U. U. = p v^ (^ H- ^) . (153.) 

1 d 



For example, if the mean flow is unidirectional, then the Reynolds stress is 



1 T Su 

T =-pu'w'=pV -^T* 

zx z oz 

where r is the stress on a surface perpendicular to Z-axis in the direction 
zx 



of X-axis. The Reynolds stress component r =-p u'w' represents the mean rate 

zx 

of turbulent transport of u' - momentixm by the w' component of the flow, where 
u', v', w', are the turbulent fluctuations of velocity components upon the mean 
flow as mentioned in equation (150). At the ocean or lake surface, the very 
top layer may be considered as Edman's "sensitive wind-driven layer." Theo- 
retically, the stress on the air-water interface should be approximately the 
same whether the measurement is made from air or from water. If the magnitude 
of surface shear stress is known and if the vertical velocity gradient can be 
measured near the interface, then the vertical eddy viscosity can be evaluated 
by using equation (I5I). 

A characteristic velocity of the turbulent fluctuating motion called the 
friction velocity is defined as 



U. = (li) V2 



., = (-^) "'" (152) 

which is a measure of the intensity of turbulence, where p is the density of 

the fluid. It can be written in terms of drag coefficient and the wind speed 
at 10 m height, 

Cio=(^)^ (155) 

Uio 



8? 

where Cio is the drag coefficient with reference to the 10 m height wind, Uio. 
The values of Cio as a function of the wind speed in m/s at 10 m level have 
been given in many empirical equations, such as (Pierson, I964): 

-1/2 
Cio = 0.009 Uio ^ (By Neimann), 

Cio = (1.00 + 0.007 Uio) 10*"^ (by Deason and Webb), 

Cio = (0.80 + O.llJi Uio) 10"^ (by Sheppard). 

In the case of a wind-driven lake current a friction velocity may be found, 

according to Deacon (I962), as 

^a 1/2 
U^ = (0.0012 —) Uio J (15^) 

^w 

where p and p are densities of air and water, respectively. 
aw 7 X i/ 

Elder and Soo (I967; have analyzed wind profile data from The University 
of Michigan research tower one mile from shore in Lake Michigan and came up 
with a best fit exponential relation of friction velocity and Uio under neutral 
conditions as shown, 

U^ = 0.066 exp (0.165 Uio) . (155) 

Following the evaluation of friction velocity, if the vertical velocity gra- 
dient is measured or estimated, the vertical magnitude of eddy viscosity can 
be calculated. 

MIXING LENGTH AND LOGARITHMIC VELOCITY DISTRIBUTION 

By analogy to the mean free path in the kinetic theory of gases, Prandtl 
introduced the "mixing length" of the fluid and von Ka'rman related it to the 
property of the mean motion. Near a solid boundary such as the lake bottom. 
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the theory of Prandtl and von Karman leads to the logarithmic law of velocity 
distribution^ 

z + z 
U-^ o / ^. 

o o 
where U^ is the friction velocity near a solid "boundary^ and z is the "rough- 
ness length," k is the von Karman constant. The logarithmic law has been used 
widely in all branches of fluid mechanics, especially in meteorology. In the 
region of neutral stability, the vertical eddy viscosity is shown (Bowden, 
V^Gh) as, 

A = k U-x- (z + z ) '^ k U^ z , (157) 

z o 00 ' V ^ I / 

since z is usually small, 
o 

Bowles et al., (1958) and Charnock (1959) used this method to investigate 

the velocity profile in the English Channel and Irish Sea and confirmed the 

approximate form of the logarithmic velocity distribution up to 2 m from the 

bottom. The value of z was from 1 to 5 mm. Lesser (l95l) measured the cur- 

o 

rent at a different level near the bottom in homogenous water and found the 
logarithmic law was valid and the roughness length for gravel-sand bottom to 
be 1.3 mm, and for mud-sand, 1.6 mm. 

Assuming the logarithmic mean velocity profile is valid in Lake Michigan 
near the bottom, then the friction velocity can be deduced from a point cur- 
rent measurement near the bottom and the vertical eddy viscosity can be calcu- 
lated. 
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EMPIRICAL FORMULA. FOR ESTIMATING VERTICAL EDDY VISCOSITY FROM WIND 

The vertical eddy viscosity can at least be approximated by the surface 
wind in a homogenous sea by (Sverdrup et al., 19^2), 

A = 1.02 W^ f or W < 6 m/s (Thorade) (158) 

A = U.3 w^ f or W > 6 m/s (Ekman) (159) 

where W is measured in m/s. 

EDDY DIFFUSIVITY BASED ON STATISTICAL THEORY OF TURBULENCE 

In a real geophysical fluid, though pure isotropic turbulence rarely 
exists, the local isotropy as introduced by the similarity theory of Kolmogoroff 
(19^1) has been proved valid. The hypothesis assumes that there exists a con- 
tinuous spectrum of fluctuations in the turbulence. The small-scale components 
of the turbulence are approximately in statistical equilibrium. They owe their 
existence to the nonlinear interchange of energy between different wave-number 
components. Apart from the effects of the variation of the viscosity and 
another parameter which is determined by the large-scale components of the 
turbulence, this equilibrium is universal. The continuous spectrum of fluctu- 
ations is considered as different scales of eddies. The intermediate scale 
eddies receive energy from large eddies (which are supposed to receive their 
energy from mean flow) and the intermediate eddies pass the energy to the small 
eddies, which finally dissipate it. Obviously, large eddies are anisotropic 
and small eddies are isotropic. Intermediate eddies can be thought of as 
qua s i - i s ot r opi c , 
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In recent years, investigations of turbulent diffusion using "marked" 
particles have become very prevalent. By introducing into the turbulent water 
such "marked" particles as Rhodamine B solution or other similar markers which 
do not affect the dynamics of the fluid and which have almost the same density 
as the fluid, the behavior of the marked particles will show only the turbulent 
characteristics of diffusion of the fluid. In this way the effect of molecular 
diffusion, which might be important as suggested by Batchelor and Townsend 
(1956), is separated out from the turbulent diffusion due to their extreme dif- 
ference in scales. In turbulent diffusion analyses both Eulerian and Lagrang- 
ian descriptions are used. If the scale of turbulence is small as compared to 
the size of the configuration of the patch of marked-particles, an Eulerian 
type of analysis is used. If the turbulent scale is much larger than the di- 
mension of the patch, a Lagrangian type of analysis is preferred. For con- 
tinuous source experiments the diffusion measurements are made relative to 
fixed positions, and the Eulerian analysis is appropriate. 

In a homogenous flow field, the turbulent dispersal of marked particles 

is the phenomenon of relative diffusion. For simplest illustration, let X,(t) 

J 

and q.(t) be the position vector and the velocity vector of the j-particle in 
J 

the marked particles patch at time t after release. Then let 

D.(t) = X.(t) - X (t) (160) 

J c 

*,(t) = q,(t) - q (t) 

J J ^ 

where D.(t) and <I).(t) are the relative displacement and velocity of the j 
particle with respect to the center of the mass of the patch, X (t), at time t. 
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If all particles are identical as assumed, then the mean relative distance is 
zero at the center of the patch. According to the central limit theorem the 
distribution of marked particles will be Gaussian after a long period of time. 
Then the average value of D.(t) is the standard deviation of the patch at time 

o 

t, and q (t) is the mean velocity, 
c 

Considering a stationary process, the auto-correlation function at two 
points can be converted to a correlation function of the time difference be- 
tween the velocities, as 






where ( ) indicates ensemble average and r = ti - t^. 

If, at t = 0, the particle is at the original position, then after suffi- 
ciently large diffusion time has elapsed, due to continuous independent move- 
ment of the marked particle, the mean variance is 



1^ r^ r^V 



D^(t) = 2 (D^ /^ /^ f(T) d tW T. (162) 



For small t, f(T) « 1. Therefore, 



D^(t) = 0^ t^ +F(oy (165) 



where D(0) is the initial separation. For large t, f(T) -> and 

A = /^ fir) d T, (164) 

will approach a limiting value which is considered as the scale of turbulence. 
Hence for large t. 



92 



D^(t) = 2 cD^ A t . (165) 

Bachelor (1953) bas proposed a development based on a similarity hypothesis of 
turbulence for discussing the relative diffusion* Its application is limited 
by the dependence of the relative diffusion on the turbulence only through the 
parameters as viscosity, Vy energy dissipation, E, diffusion interval time, t, 
and the initial separation, D(0). Dimensional arguments, then show that 
(Okubo, 1962) 



d D (t) , / X x-> /q 

= pi(E D(0)r/^ t , for small t. 



(166) 



d D (t) p 

— — = P2 E t , for intermediate t (167) 



and 



— -^ — = P3 ^^ A asymptotically, (168) 

where pi, p^, P3 are constants. 

By analogy to the development of the molecular diffusion when the scale 
of dispersion becomes large as compared with the largest eddy present, the 
equivalent eddy diffusivity may be shown as 



1 d D^(t) ^ .. , 



or 



1 D^(t) , 
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RICHARDSON'S NEIGHBOR SEPARATION 

Richardson (1926) and Richardson and Stommel (19^8) have considered 
"neighbor separation." Instead of the diffusivity K, they use F(£), the "dif- 
fusivity for neighbors." The ordinary Pick's diffusion equation is replaced by 

where q(i) is the neighbor concentration function, which is the number of 
neighbors per unit projected length It. From observed atmospheric data, Richard- 
son deduced the well known k/'^ power law of relative diffusion, 

F(i) = C ^V3 ^ (1^1) 

where C is a constant. 

In the ocean as well as in the Great Lakes, a thermocline usually exists 
and migrates up and down seasonally. The induced vertical gradient in density, 
mostly in stable situations, causes a reduction of vertical turbulent velocity 
by converting a portion of kinetic energy to increase the potential energy. 
Thus, in a stable vertical density gradient, the intensity of turbulence is 
reduced, hence the eddy viscosity, A , and eddy diffusivity, K , are also re- 
duced. The latter will decrease more than the former (Hill, I962). In the 
sea as well as in the Great Lakes, the horizontal components of turbulence are 
usually much greater than the vertical component, not only because of the vast 
difference in dimensions but because of the effects of the vertical density 
gradients. Hence the horizontal eddy viscosity and the eddy diffusivity are 
expected to be much greater, though the horizontal and vertical mixing trans- 
port of the fluid property may be still of the same order due to much larger 
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vertical gradient of the fluid property (Defant, l-96l), Csanady (I963) has 
shown that the vertical eddy diffusivity of the Great Lakes is only one or two 
orders of magnitude greater than molecular viscosity (about 0.1 - O.5 cm^/sec) 
while horizontal diffusivity is 10^ to 10^ times greater. Moreover, the lon- 
gitudinal horizontal diffusivity is 4 to 5 times greater than the lateral dif- 
fusivity. Bowden (1965) has shown K (longitudinal diffusivity) is about 10 
times greater than K (lateral diffusivity), and the mixing produced by a 
shearing current becomes important if a moderate degree of stability is present. 
In the Great Lakes, a shearing current has been observed to increase the effec- 
tive horizontal diffusivity by a factor of 5 or so (Csanady, I966). Foxworthy 
and Barsom (I967) have come to the same conclusion as Csanady that increasing 
stability results in a lower vertical diffusivity, K , causing a more complex 

current pattern. This in turn increases the horizontal viscosity, K , and the 

«y 

net result is that the product of K K is nearly constant. 

z y 

The Magnitudes 

Eddy viscosity and eddy diffusivity are essential parameters for the 
theoretical thermal current model study of Lake Michigan. Therefore, a series 
of experiments were carried out to evaluate the representative order-of-magni- 
tude values for the eddy viscosity and eddy diffusivity to be used in the mod- 
eling investigation. 

AN ESTIMATE OF EDDY VISCOSITY BY REYNOLDS STRESS AND VERTICAL VELOCITY GRADIENTS 
In August 1967 a pair of current meter buoy stations were set 8 miles off- 
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shore from Sheboygan^ Wisconsin. During the period of the current measurement^ 
the research vessel R/V INIAND SEAS was operating around the stations (during 
the daytime) whenever weather permitted. Winds ^ as well as other weather data^ 
were automatically recorded on a digital meteorological tape recorder. The 
anemometer and vane that fed into the recorder were on the mast of the INIAM) 
SEAS 1^1- m above the water surface. The Captain's log and hourly meteorological 
observations recorded wind from a lower vane and anemometer at about 11 m above 
the surface. Using the analyzed recorded mean wind versus height, plotted on 
a log scale with the mean velocity zero at about 1 cm height, (Floyd C. Elder, 
personal communication) the 10-m level wind was interpolated with reference to 
the Captain's log. The 10-m daily winds are shown in Table k. The vertical 
velocity gradients in the direction of the surface current were deduced from 
two current meters separated vertically by 10 m. According to equation (15U), 
the Reynolds stress is 

T = p ul = 0.012 p U?o . 
t w a 

Then using p = 1 g/cm^ and p = 1.23 x 10'^ g/cm^ the eddy viscosity can be 
w a 

calculated. The data, wind, wave conditions, absolute vertical velocity gra- 
dient, shear stress, and calculated eddy viscosity are shown in Table h. The 
values of Reynolds stresses obtained here by using Deacon's formula are in good 
agreement with the mean value of stress obtained by Elder and Soo (1967) at 
The University of Michigan Research Tower on the east side of the lake. 
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AN ESTIM/VTE OF THE UPPER LIMIT VALUE OF EDDY VISCOSITY AND EDDY DIFFUSIVITY BY 
STATISTICAL THEORY 



It is quite well known (e.g.^ Ayers, et al,, 1958) that a large clockwise 
eddy exists in southern Lake Michigan between Grand Haven and Michigan City. 
This large eddy has the length of scale of a little less than half the dimen- 
sion of the width of the lake. Bowden (V^Gh) has shown the effective eddy 
viscosity is proportional to the 4/5 power of the length scale L, of the par- 
ticular eddy under consideration, similar to the Richardson's formula equation 
(171). In Lake Michigan^ let this largest eddy have an average length scale 
of L = 40 km or L = 4 X 10^ cm. Then, according to Richardson's formula, the 
eddy visQogity will be of the order lO'^ cm^/sec which is the upper limit of 
the numerical values of th^ eddy viscosity and eddy diffusivity in Lake Mich- 
igan. 



AN ESTIMATION OF VERTICAL EDDY VISCOSIT^f FROM CURRENT MEASUREMENT NEAR THE 
BOTTOM 



From 11 May to 11 October 19^1, a tripod-supported pendulum current meter 
was installed under the direction of Dr. Ayers, about 25O m off shore, 7-1/2 
miles south of Benton Harbor at a depth of 5 m. The pendulum was suspended 
only 20.5 cm above the sandy bottom. The monthly mean current is shown in 
Table 5. Using Lesser 's (l95l) roughness length for sand bottom, z = .15 cm, 
the friction velocity near the bottom can be determined through equation (I56), 

z 

U^ = k U In ( 2—) 

z + z 
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where U is the mean velocity, z is the height from the bottom during measure- 
ment, and k = 0.4l. The monthly mean friction velocity for z = 20.0 cm is 
o 

shown in Table 5. The neutral stability condition generally persists near the 
bottom. Then the vertical eddy viscosity can be estimated by equation (15T) 
as shown in Table 5 • 



TABLE 5 

Estimated Vertical Eddy Viscosity from 
Current Measurement Near Bottom 





Mean 


Friction 




Vertical eddy 


Month 


current 


velocity 




viscosity 




(cm/sec) 


(cm/sec) 


(Az 


; = k u*z cm^/sec) 


May 


59.5 


k.Gk 




Uo.o 


June 


31.7 


2M 




20.7 


July 


kl.2 


5.68 




30.9 


August 


58.5 


i^.56 




38.3 


September 


56.1 


1^.38 




36.0 


October 


69.8 


5.44 




h5.1 


Average 


; value 






55 cm /sec 



AN APPROXIMATION OF THE VERTICAL EDDY VISCOSITY BY SURFACE WIND 

Wind data from both the two-current-meter-buoy operation off Sheboygan 
and the wind and current station outside of the Benton Harbor beach are used 
for approximating the vertical eddy viscosity in Lake Michigan, In the two- 
current -meter-buoy operation off Sheboygan, wind data were measured aboard the 
R/V inland seas. At the current meter station outside of the Benton Harbor 
beach, the wind was measured Just above the water. From Sverdrup et al., 
(19^2), depending on whether the magnitude of wind was greater or less than 
6 m/s, the vertical eddy viscosity was approximated by equations (158), (l59) 
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as shown in Table 6. Notice that the vertical eddy viscosity is smaller when 
the wind is averaged for a longer period of time. 



TABI.E 6 



Estimated Vertical Eddy Viscosity from Wind Speed 



Time 



Wind (m/s) 



Vertical eddy 
viscosity, A cm^/sec 



Remarks 



IT Aug. 


1967 


9.5 


18 Aug. 


1967 


9.8 


20 Aug. 


1967 


5.7 


21 Aug, 


1967 


3.h 


23 Aug. 


1967 


6.2 


2k Aug. 


1967 


k.6 



Average value 



381 

413 
140 

ko 

165 

100 



Wind averaged 
over 6 hours 



206 cm^/sec 



May 1967 



5.5 



June 1967 


2.5 


July 1967 


3-1 


Aug. 1967 


3-2 


Sept. 1967 


2.8 


Oct. 1967 


i^.l 


Nov. 1967 


k.k 



Average value 



16 
30 

33 
22 

70 

87 



Monthly mean 
wind 



43 cm /sec 



MEASUREMEM'S OF THE EDDY DIFFUSIVITY BY THE DISPERSION OF MARKED PARTICLES 

In a series of experiments conducted in Lake Michigan, the fluorescent 
dye rhodamine B was used as a tracer for diffusivity measurements. Rhodamine 
B has a peak fluorescence wavelength of 58O m|a. In May I968, a fluorometer 
with only one sample cell was used. At this season the phytoplankton was 
abundant, and since they have an ultraviolet fluorescence frequency Just a 
little lower than that of rhodamine B, the fluorescence of the peak concen- 
tration of rhodamine B dye was not distinguishable from the background after 
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an hour or so. In all subsequent experiments, a comparative measuring fluoro- 
meter with two cells was used. The rhodamine B mixed with natural lake water 
flowed through the dye sample cell, while the standard cell contained natural 
lake water. The relative difference in fluorescence between two cells was 
measured to obtain reliable data in the possible presence of any intensity 
variation of the ultraviolet light lamp, systematic changes in the character- 
istics of photocell, and differences in the amount of naturally occurring flu- 
orescent substances in the water (Noble and Ayers, 1961). Both instantaneous 
source and continuous source measurements were made in I968. 

Instantaneous Sources 

Three runs of instantaneous sources of marked particles for diffusivity 
studies were conducted. These were 1 and 8 miles off Grand Haven in May, and 
5 miles off Milwaukee in June. Rhodamine B dye was diluted with methanol to 
have a specific gravity very close to the lake water. Generally, one part of 
stock solution of rhodamine B in kO^o acetic acid solution was mixed with 5-1/2 
parts of methanol. The dye mixture was gently poured out on the surface lake 
water at the upwind side of the boat in order to allow the boat to drift away 
with least interference to the dye patch. The first and second runs were in- 
vestigated aboard the R/V MYSIS and the third run aboard the R/V INLAND SEAS. 
During the first 20 min after the source had been instantaneously released, 
the diameters of the expanding dye patch were estimated by eye with the ship^s 
length at reference. At 20-min intervals, the boat was used to coast through 
the dye patch longitudinally and laterally at the lowest possible speed. For 
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each transect, fixes were obtained from shore objects when entering and leaving 
the patch, along with the time required to traverse through the patch. Running 
back and forth across the dye patch, the concentration of the dye patch was 
continuously recorded from the fluorometer on an AZAR recorder. The concentra- 
tion distribution of a typical pass through the patch is shown in Figure 9. 




IN 



2^ 3' 4' 

Distance (in unit scale) 



6'OUT 



Figure 9« A typical distribution of dye concentration 
across a dye patch. 



Most generally, the dye concentration distribution showed a form of Gaussian 
distribution. The diameter of the patch, longitudinally or laterally, is prop- 
erly approximated as k times the standard deviation of the patch from its mean 
position. Batchelor (1955) shows that (assuming the velocity field is uniform, 
homogeneous, and steady) at the initial stage of a diffusion process, the ef- 
fective diffusivity (which is proportional to the rate of change of the square 
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of standard deviation with time) is proportional to the first power of the 
elapsed time; in the intermediate phase, it grows as a 4/3 power of the patch 
size; and at the final stage, as the standard deviation grows as 1/2 power of 
the elapsed time, the diffusivity should symptotically approach a constant. 
The eddy diffusivities have been evaluated according to equations (l69) or 
(170). The growth of dye patch versus time is shown in Figures 10 and 11. 
General information and values of diffusivities are shown in Table 7. 



TABLE 7 

General Information and Effective Eddy Diffusivity 
of Dye Patch Studies 





Run 1 


Run 2 




Run 5 


Date 


1 May 1968 


2 May I968 




27 June 1968 


Location 


1 mile off shore 


8 miles off 


shore 


h miles off shore 




near Grand Haven 


near Grand Haven 


near Milwaukee 


Depth 


20 ft 


If-O ft 




65 ft 


Duration of run 


1130 - 1210 EST 


1300 - 1450 


EST 


0910 - 1240 EST 


Current velocity 










(cm/sec) 


11.6 


5.5 




7.6 


Wind 










Direction 


m 


SW 




ssw 


Speed (kt) 


13 


15-20 




5 


Lake Condition 


2-3 ft waves 


1 ft wave 




2 ft waves 


Effective eddy dif- 










fusivity (cm^/sec) 










Lateral 


Not available 


1140 




7k3 


Longitudinal 


Not available 


4731 




4019 


Remark 


Very poor visi- 
bility 
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Continuous Source 



Four runs of continuous (in time) point sources of marked particles were 
conducted near the meteorological research tower 1 mile off the east shore in 
Lake Michigan • Ten gallons of rhodamine B and methanol mixture were contained 
in two 5-gSil carboys connected together by plastic tubing as shown in Figure 12, 



Air Intake 




Screw Clamp 



Figure 12. Continuous source apparatus. 



The continuous point source was released from the apparatus through a rub- 
ber tube which had its outlet about 20 cm above the water surface. The rate 
of discharge of the dye was controlled by a screw clamp on the tubing, and was 
usually allowed to flow at a rate of J gal/hr, or approximately k gm/sec. In a 
homogeneous fluid with a steady uniform current, the continuous dye source 
produced a slender plume extending downstream. 



io6 

The direction of the dye plume was reported by the observer using a theod- 
olite mounted on the meteorological research tower 57 feet above the water sur- 
face. Buoys were dropped as markers on the edges of the dye plume at right 
angles to the direction of flow. The approximate distance of the marker from 
the source depended on the mean velocity of the current which was estimated 
from drogues set earlier. 

The R/V miosis crossed the dye plume back and forth at 10-min intervals at 
the lowest speed she could make (about ^.5 mph). As soon as the dye plume 
passed through the markers, samples of water were taken through a constant-flow 
fluorometer that fed into an AZAR recorder. Communication between the tower 
and the boat was maintained through a VHF transceiver, "Whenever the boat en- 
tered and left the dye plume a "MARK" signal was transmitted to the theodolite 
observer. The observer obtained the azimuths of both edges of the plume and 
the angular elevation from the tower to the boat. 

The horizontal distance between the boat and the tower was obtained from 
the elevation angle as well as by fixes from land objects. The angle subtended 
by the edges of the plume is the difference of the two azimuths. The arc 
length of the plume at this distance from the tower was arbitrarily taken as 
four times the lateral standard deviation unit. At each return crossing, the 
vertical diffusion was measured by lowering the fluorometer intake at half- 
meter intervals until there was no indication of the dye. Half of the maximum 
vertical displacement was taken as the vertical standard deviation unit. 

In the last few minutes before the exhaustion of the point source, con- 
tinuous crossings were conducted at various distances along the plume from the 
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tower to approximately 1-1/2 miles downstream, These crossings measured the 
asymptotically steady diffusion of the plume. 

These continuous point source diffusion studies were carried out on JO 
July and 1-2 August 1967- On 50 July, sampling was discontinued about 2 hr 
after the experiment was initiated, because of high seas. The dye plume was 
continuously sampled for 5 hr on 1 August, and during the sampling period the 
current shifted its direction of flow about 15"". Two sampling runs were con- 
ducted on 2 August under a very steady current condition. The general infor- 
mation for each run and the values of diffusivity obtained are shown in Table 
8. The lateral growth of the dy^ plume versus tim^ at a constant distance 
from source is shown in Figure I5 . The lateral growth versus distance under 
asymptotically steady condition is shown in Figure l4. 

A MEASUREMENT OF THE HORIZONTAL DIFFUSION BY DROGUE STUDIES 

A drogue is made of h pieces of light sheet metal (e.g., aluminum alloy) 
approximately 1x2 m in dimension, linked togeijiher in the shape of a cross 
and suspended by chains in the water by a float which supports a radar reflec- 
tion above the water surface (similar to Farlow, I965 ) . It is frequently used 
for measuring the mean current of the flowing layer at the drogue depth. The 
equivalent radius of the drogue is about 1 m and, based on its geometry and 
the corresponding Reynolds numbers (lO^ to 10^ in the Great Lakes), its time 
constant is estimated to be 100 sec (Okubo and Farlow, I967). Though the 
small-scale turbulence may be averaged out, the drogue is adequate for measure- 
ment of large-scale turbulent diffusion for eddies larger than 10 m in scale 
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which represent the energy- containing eddy scales in the Great Lakes (Csanady, 
1963; Okubo and Far low I967). 

On 27 June I968 a hexagonal pattern of surface drogues, with the central 
one very close to a fixed buoy, was set about 5 miles off shore near Milwaukee 
in Lake Michigan. The original distance between each pair of drogues was 1/2 
mile- After setting the pattern, the R/V INLA.NI) SEAS continuously fixed each 
drogue for the following 6-1/2 hr. Similar to the result mentioned by Okubo 
and Farlow (1967), the growth of the pattern size with time was not clearly 
noticed* The diffusion constant was deduced from the second moments about the 
mean of the hexagonal pattern. At any particular time, the instantaneous po- 
sitions of all drogues were interpolated from their respective previous fixes. 
The general information on the experimental run and the value of eddy diffu- 
sivity are shown in Table 9, 

TABLE 9 

General Information and Eddy Diffusivity 
of Drogue Study in Lake Michigan 



Date 27 June I968 

Location 5 miles off Milwaukee 

Depth of water 65 ft 

Depth of drogues 5 ft 

Duration of run 1207 - 183 EST 

Current velocity 7,6 cm/sec 

Wind 

Direction N 

Speed 6 kts 

Lake Condition 1/2 - 1 ft waves 

Mean initial separation 2378 ft 

Final deviation from mean 5^00 ft 

Effective eddy diffusivity 3.4 x 10^ cm^/sec 
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Comparisons of Values with Other Studies 

In Lake Michigan, there is little data about eddy viscosity or eddy dif- 
fusivity. Some figures from the ocean or other lakes are available as refer- 
ence. In modeling the Gulf Stream, Munk (1950) estimated the lateral eddy 
viscosity as 5 x lo"^ cm^/sec. Stommel (1955) estimated the eddy viscosity of 
the Gulf Stream near the Florida Straits to be 2 x 10^ cm^/sec. Defant (196I) 
determined the exchange coefficient of lateral large-scale turbulence in ocean 
phenomena as 5 x 10^ cm^/sec. 

In general, the vertical turbulence scale is much smaller in dimensions 
than that of the horizontal. Therefore, the magnitude of vertical eddy vis- 
cosity is much smaller. According to Defant (1961), Sverdrup obtained a value 
for the vertical eddy viscosity of 75-260 cm^/sec in the North Siberian Shelf 
region; Fjeldstad estimated the vertical eddy viscosity in the Caspian Sea as 
0-22i^- cm^/sec, and Suda measured the values 68O-75OO cm^/sec for the Kuroshio 
and 15 0-146 cm^/sec in the Japan Sea. 

For vertical eddy diffusivity, some typical values have been summarized 
by Defant (1961) as: Mediterranean, k2 cm^/sec; California current, JO-UO 
cm^/sec; Caspian Sea, 2-l6 cm^/sec; Equatorial Atlantic Ocean, 52O cm^/sec. 
Csanady (1963, 1964, I966) has measured the eddy diffusivity in Lake Huron and 
Lake Erie, and found horizontal diffusivity of the order of 5OO-25OO cm^/sec 
in Lake Huron and IOOO-5OOO cm^/sec in Lake Erie. These values are much greater 
than vertical eddy diffusivity, which is only of the order 0.1-10 cm^/sec. 
Okubo and Farlow (1967) measured the horizontal diffusivity for large eddies 
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in Lake Michigan and Lake Erie using drogues and found the eddy diffusivity to 
he 3 X 10^ to 6 X 10^ cm^/sec. 

According to our estimation^ the vertical eddy viscosity is on the order of 
1.2 X 10 to 2 X 10^ cm^/sec in Lake Michigan. The lateral eddy diffusivity we 
found using marked particles is 750-1200 cm^/sec for small eddies and ^ .k x 10^ 
cm^/sec for larger eddies. The effective longitudinal eddy diffusivity for 
small eddies is in the order of ^ x 10"^ cm^/sec. All the data we used appear 
quite consistent and in good agreement with that of others. 

As Csanady (1966) has pointed out^ shearing flow and unsteady currents or 
any kind of complex circulation patterns accelerate the diffusion very repidly. 
Therefore^ the above estimations may still be in the lower range of eddy vis- 
cosity or diffusivity in Lake Michigan. 

Reasonable Values for Model Study 

In modeling a geofluid problem, the eddy viscosity and the eddy diffusivity 
are of critical importance in order to predict the natural current or wave 
phenomena. In Lake Michigan, though the mean current velocity is small in 
general, the flow field is turbulent in nature. In solving or explaining the 
flow pattern of the mean lake current, it is possible to use the governing 
equations of laminar flow with eddy viscosity and eddy diffusivity in place of 
the molecular viscosity and the molecular diffusivity, A series of experiments 
have been conducted in order to estimate these two values. The data presented 
are quite consistent and in good agreement with the data reported by other in- 
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vestigators. The vertical eddy viscosity in Lake Michigan is in the range of 
1 to 10 with a mean 50 cm^/sec. The horizontal viscosity is in the range of 
10^ to 10^ with a mean value of 10^ cm^/sec. The eddy diffusivity may reach 
the same magnitude as the viscosity but it is in general smaller. A typical 
mean value for vertical eddy diffusivity is 10 cm^/sec^ and for the horizontal 
eddy diffusivity 10^ cm^/sec. Though the few experiments are far from enough 
to state the characteristic values of turbulent phenomena, the data presented 
together with values reported by others can be used as good reference values 
for Lake Michigan. 



VI. REPRESENTATIVE VALUES OF PARAMETERS 
IN THE LAKE MICHIGAN MODEL 



General Physical Data 



General features of Lake Michigan have been discussed in detail in Chapter 
I and II. 

All physical parameters used for our theoretical model study in similitude 
to the real phenomena in Lake Michigan are summarized in Table 10. Note that 
all notations have been specified in Table 2. The magnitudes of e, £i, e^ 
should also be noted. 

TABLE 10 

General Physical Data Used in 
the Lake Michigan Model 



Parameter 
Notation 


Value 


Unit 


L 


1,22 X 10^ 


cm 


D 


1.22 X 10* 


cm 


2Q 


10-* 


rad/sec 


V 
V 


50 

5 X 10* 


cm/sec 
cm/sec 


10^ 




e 


2 X 10"® 




7 


10"^ 




a 


1 




K 


11 




El 


k X 10"^ 




£2 


k.k X 10"^ 




s 


980 


cm/sec 
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Parameters as Functions of Seasonal Temperatures 

The horizontal temperature difference AT as well as the thermal Rossby 
number^ R, are functions of seasonal temperatures. The reference velocity V, 
satisfying the geostrophic thermal wind relationship, varies as the seasonal 
temperature changes . 

GENERAL THERML CURRENT 

In the type-A or type-B circulations the equation of state is approximated 
as equation (7). The coefficient of thermal expansion of freshwater density, 
Qy is a special function of local temperature. Since the whole lake is above 
if°C, during the summer-heating period in the type-A circulation and autumn- 
cooling period in the type-B circulation, a is positive in equation (7), that 
is, as the local temperature T' increases away from the temperature of maximum 
density, the density of water decreases. The horizontal temperature difference, 
AT, ranges from 0° to 10°C. Hence, the coefficient of thermal expansion of 
fresh water near 11°C is taken as a representative value. All values of special 
parameters connected with seasonal temperature changes used for general thermal 
current are listed in Table 11. 

During the winter-cooling and spring-heating periods the temperature in 
the lake is all below 4° C. The horizontal temperature differences during these 
periods are smaller due to the occurrence of icing. The coefficient of thermal 
expansion, a, is negative in the equation (7) as the local temperature, T', is 
less than T , but the density of water is still decreasing as T* deviates away 
from k^'C. The approximate general thermal- cur rent circulations during these 
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TABLE 11 



Values of Parameters Connected With Temperature 
for the General 'niermal Currents 

Value 

Spring-heating Winter-cooling 
Jrarametier 

. . . and autumn- and autumn- Unit 
notation 

cooling periods heating periods 



a 


2 X 10""* 


2 X 


10" 


■5 


°c-^ 


AT 


10 


k 






°c 


V 


20 


1 






cm/sec 


R 


1.65 X 10"^ 


10"^ 









periods are deduced from the circulation pattern during the summer -heating and 
autumn-cooling periods . It can easily be checked out either physically or 
mathematically that the winter-cooling period and the summer-heating period 
possess the same circulation pattern due to the combined effects of the special 
character of the a above or below the temperature of maximum fresh water den- 
sity and the opposite imposed horizontal temperature gradients. For similar 
reason, the period of autumn-cooling and the period of spring-heating have the 
same circulation pattern. Therefore, with the circulation pattern of the sum- 
mer-heating period as reference, we can still take a as positive, but treat 
the horizontal temperature difference in the negative sense during winter- 
cooling and spring-heating periods. The coefficient of thermal expansion at 
3°C is used as representative for these two periods. 

All values of special parameters connected with seasonal temperature 
changes used for general thermal current are listed in Table 11. The magnitude 
of R should be noted. 
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THERML BAR CURREINT 

With the same reasoning, either the spring thermal bar or the autumn ther- 
mal bar periods have the same circulation pattern except the latter is much 
weaker. The equation of state during thermal bar periods is approximated by 
equation (89). All values of parameters connected with temperature during 
thermal bar periods are listed in Table 12. The magnitude of R should be noted. 

TABLE 12 

Values of Parameters Connected 

with Temperature for 

the Thermal Bar Currents 



Parameter 
Notation 


Values 


Unit 


A 


4 X 


10-3 


°c-^ 


AT 


8 




°C 


V 


12 




cm/sec 


R 


10"2 







VII. THERML CUEEENT CIRCUIATION PATTEMS AND TEMPERATUEE 

DISTRIBUTIONS FROM SOLUTIONS OF THE BAROCLINIC 

DYNAMIC MODEL OF LAKE MICHIGAN 



In our simple categorization for thermal current circulation as functions 
of seasonal temperature differences as introduced in Chapter I^ the type-A 
circulation pattern induced during the winter-cooling period and the summer- 
heating period^ and the type-B circulation pattern induced during the spring- 
heating period and the autumn-cooling period, belong to the general thermal 
currents. The thermal bar circulation pattern induced during both spring and 
autumin thermal bar periods is the type-C circulation in our terms. The surface 
temperature profiles in different periods have been mentioned in the above 
chapters. Note also that the strong thermal stratification period is not in- 
cluded in our- quasi-homogeneous model study. 

The general thermal current circulation pattern and temperature distribu- 
tion have been mathematically investigated in detail in Chapter III, using the 
summer-heating features as representative. The special thermal circulation 
pattern during the thermal bar period has been demonstrated theoretically 
in Chapter IV. Following the double perturbation procedures, with regular 
perturbation expansion in powers of the thermal Rossby number R, and then sin- 
gular perturbation analysis with respect to the Ekman number e for each power 
of R, analytic approximate solutions of the whole flow field have been obtained 
except in the small corners of the cross-section. Representative values of 
parameters closely in similitude to the real field phenomena in Lake Michigan 
have been put into the solutions obtained in Chapters III and IV by computer. 

119 



120 
The results shall be discussed in more detail below. 

Temperature Distribution 

GENERAL TEERmL CURRENT PERIODS 

The zeroth R-order temperature is shown in equation (37) which has no 
boundary layer effect, and the first R-order temperature is shown in equation 
(71) with boundary layer modifications in each individual boundary region as 
shown in equations (75 )> (T^)^ (75 )> and (76). According to the double per- 
turbation method (equations (I7) and (5I)), the total temperature, is 



T = T^°) + Ral T^'^ + ei^/^(iT^') + ^T^'^ 



+ £2'-/^(3T^^^ + 4T^''^)] + O(cr^R^). (I72) 

Remember the boundary solutions are only of local importance at their respec- 
tive boundaries . We notice from the above equation that the boundary layer 
effect of temperature is not very important. 

Using representative values of parameters given in Chapter VI, the non- 
dimensional temperature distribution is shown in Figure I5 . During the summer- 
heating period, the surface temperature is concave cosine profile and the hor- 
izontal temperature is much higher on the edges than in the middle portion of 
the lake. We consider this kind of temperature difference with higher tempera- 
ture near the shore as a positive gradient. With data given in Table 11, each 
isotherm in Figure I5 is about 1°C in difference. Though the middle portion 
of the lake may still be at ^°C, the two edges of water near shore may have 
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been warmed up to l4°C. 

During the autumn-cooling period, since the atmospheric temperature is 

colder, the surface temperature of the lake has a cosine profile with the two 

edges at nearly the reference temperature, T . The horizontal temperature dif- 

o 

ference is in the negative sense because the edges of water are colder. The 
temperature profile has exactly the same form in isotherm distribution, but 
with a reversed direction of temperature difference as compared with that of 
the summer-heating period, that is, the temperature gradient is in the negative 
sense with maximum temperature at the middle during this period. 

For the other two periods of time under general thermal currents the tem- 
perature isotherms have the same form except that the magnitudes of temperature 
differences are smaller. During the winter-cooling period, the maximum surface 
temperature is nearly the temperature of maximum density in the middle of the 
lake while the edges may well be frozen. The temperature distribution is of 
exactly the same form as that of the autumn-cooling period except the tempera- 
ture difference is only k°C. The spring-heating period is exactly the same 
as the summer-heating period except the horizontal difference is smaller. 

THERJVIAL MR PERIODS 

The zeroth R-order temperature of the thermal bar period is shown in equa- 
tion (109) and its first R-order overall part of temperature is shown in equa- 
tion (155) with 0{e'^/^) boundary parts respectively shown in equations (125), 
(126), (127), and (128). The total nondimensional temperature, with values 
of parameters shown in Tables 10 and 12, is shown in Figure I6. During the 
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spring thermal bar period, the surface temperature in the middle of the lake 
is at 2°C to 3°C but the edges near shore have already warmedup to a tempera- 
ture far above k'^C. The thermal bar exists near the zone of temperature of 
maximum density. The maximum horizontal temperature difference between the 
thermal bar and the shallowest edges may be as large as 10°C. 

The autumn thermal bar has the reverse distribution with the temperature 
in the middle of the lake slightly above 4°C and two edges of the lake much 
colder at nearly 2°C. The temperature isotherms have the same form but temper- 
atures decrease shoreward and the horizontal temperature difference is only 2 

to yc. 

Velocity Field 

GENERAL THERMAL CURREOTS 

Meridional Velocity 

Under the basic assumption that the geostrophic and thermal gradient 
forces are in equilibrium in the flow field, the meridional velocity dominates 
the circulation pattern. The zeroth R-order meridional velocities of the type-A 
circulation pattern, the circulation pattern occurring during the winter-cooling 
and summer-heating periods, are shown in equation (58) for the interior region 
and equations {k9) , (53)> (5?)^ (59) for the respective boundary regions. The 
first R-order meridional velocities are shown in equation (78) and equations 
(80), (82), (84), (88). The total meridional velocity, in form similar to 
equation (I72) for temperature, is then. 



V = V 

o 
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o 

+ O(R^). (173) 

The nondimensional total meridional velocity of the type-A circulation has been 
plotted as shown in Figure I7. In the upper half of the fluid in the Lake 
Michigan model, there is a northbound flow on the eastern part of the lake and 
a southbound flow on the western part of th,e lake. In the middle portion of 
the lake, the meridional velocity is very weak. It beomces stronger away from 
the center and reaches a maximum at places around pne-fourth of the lake width 
from shores. Then it decreases toward the shore but becomes strong again within 
the eastern and western boundary regions due to the reinforcement of the bound- 
ary layer currents near "the side slant boundaries. In general, the magnitude 
of meridional velocity first increases then decreases as depth increases. The 
meridional velocity has only a slight change in the upper Ekman layer due to 
the effect of the boundary l^yer current. The maximum meridional velocity in 
the vertical section occurs at few meters below the free surface and its mag- 
nitude may reach as high as I7 cm/sec during the summer-heating period. In the 
deeper part of the lake, the meridional velocity is much smaller. Below two- 
thirds of the lake depth, the current direction may reverse with southbound 
flow in the east and northbound flow near the west shore. 

Meridional current during the winter-cooling period possesses the same 
circulation pattern as shown in Figure I7 except the magnitude is smaller. 
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The maximum magnitude may reach 1 cm/sec in winter. 

The type-B circulation has exactly the reverse circulation pattern from 
that of the type-A circulation. 

Zonal and Vertical Velocities 

During the summer-heating period, the stream function for zonal and ver- 
tical velocity components is 0(£). The zeroth R-order stream function in the 
interior region is independent of the vertical coordinate as shown in equation 
(39) • According to the velocity definition in equation (ij), the zonal veloc- 
ity component is 

and the vertical component is 

Therefore, there is no zeroth R-order zonal velocity in the interior. The 
zeroth R-order boundary stream functions are shown in equations (50), (55), 
(58), (60). The first R-order stream functions are shown in equation (79) 
for the interior region and in equations (8I), (85), (85), (88) for the respec- 
tive boundary regions. According to equation (l7^), the total zonal velocity 
can be written as: 

u = R£i -^^ (cosh az - G) sin knx + £iV2{E/2 e"^^^^ sin -^ - — e"*^^^^ 
2 \[yE sf2 42. 

( li ^ . liM . o ^ Rjt^(l - G) -^//2 . C2 . » . 1/2 
• (cos -^ + sm -^)] li sm 2Ttx + ^ e ^' sin -^ sin 4jtx) + eV^ 

n/2 ^2 n/2 7H ^2 
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• {^ [sin 2^j{l . r^) e'^^^^ sin ^ - sin 2.1X7(1 - nJe^^^^ sin ^ ] 

nTk n/'2 ^f2 

+ R --^^ [(cosh aril - G) sin hnky{l - T)i)e^^^ sin — + (cosh ar^ - G) 

Nr2K7H ^2 

' sin 43t\7(l - Ti2)e"^^^^ sin -^ ]} + 0{R^e^/^). (I76) 

n/^2 

The total vertical velocity, referring to equation (I75), can be written as 

w = Si2Tr^{[-^ (cosh az - G) cos Uitx - cos 2atx] + [e"^^^^ cos -^ + e""^^^ 
'/7H \/2 

/ ^2 ^ . C2\i ^ ^ Rjt(l - G) -t2//^2 , ^2 . Csn . . 

• (cos -^ + sm -^jjcos 27TX + e '^ (cos -^ + sin -^) cos hjtx] 

4 2 42 E\fj 42 42 

- £1 — {[e^^^ sin — sin 2.1X7 (l - t)i) + sin 2.1X7(1 - r^) e ^^^ 
4kj 42 

• sin -^ ] - — ^ [(cosh arji - G) e ^ sin ^ sin ^.1X7(1 - tji) - (cosh aT]^ - G) 

\f2 -^ 42 

• sin 4TrX7(l - rfe) e'^^^^ sin -^]) + 0{R^e^/^). (I77) 

\f2 

With data from Tables 10 and 11, equations (I76) and (I77) have been calculated 
by computer and the results are shown in Figure I8. There are two relatively 
extensive cells which meet at the middle of the lake. The convergence of cells 
results in a broad sinking motion in the central portion of the lake and up- 
wellings at both east and west sides of the lake. The velocity components in 
the interior region are mostly vertical with only an 0(Re) part zonally. In 
the upper and lower Ekman layers, zonal components are dominant but they are 
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entirely within the surface and bottom boundary regions. On the slant side 
boundaries, the zonal and vertical velocities are much stronger than in any 
other place in the flow field. Generally, the magnitude of the mean zonal 
velocity is less than 1 cm/sec with a maxiraum near this value at the side slant 
boundaries. The mean vertical velocity is, in general, in the range of 10""^ 
to 10"^ cm/sec and may reach a maximum mean value of 2 x 10""^ cm/sec near the 
side boundaries . 

As mentioned before, the zonal and vertical velocity components of the 
winter-cooling period have the same circulation pattern as shown in Figure 18 
but the magnitude is about a factor of ten smaller. 

The type-B circulation is just the reverse of the pattern as shown in 
Figure 18. 

THERMA.L BAR CURRENT 

Meridional Velocity 

The zeroth R-order meridional velocities of the type-C circulation are 
shown in equation (112) as the interior solution and equations (ll5), (ll?), 
(119), (123) as respective boundary layer solutions. The first R-order merid- 
ional velocities are shown in equation (157) and equations (l^+O), (m2), (l44), 
(li|8). The result has been plotted in Figure I9. As we expected, the circu- 
lation pattern is more complicated. There are broad shearing flows along each 
thermal bar on each side of the lake near shore. The meridional flows in the 
main body of the lake between the two thermal bars, southward flow near the 
eastern thermal bar and northward flow near the western bar, form a clockwise 
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couple while the middle portion of the lake is a nodal point. The meridional 
velocity at the thermal bar is very weak. Meridional flows between each ther- 
mal bar and its nearest shore are northbound on the eastern side and southbound 
on the western side of the lake. The magnitude of the meridional velocity, in 
general, decreases as the depth increases. 

On the upper Ekman layer, the magnitude of the meridional velocity is al- 
most uniform. The meridional flow is very weak near the bottom. The direction 
of flow is the same through the whole depth. Probably maximum meridional ve- 
locity in spring time is about 12 cm/sec. 

The spring thermal bar current is stronger than the autumn thermal bar 
current which has a probable maximum meridional velocity of less than 1 cm/sec. 

Zonal and Vertical Velocities 

According to the stream functions given in equation (llj), (ll6), (ll8), 
(120), (124) for the zeroth R-order and in equations (139), (l4l), (Uj), 
(ll+6), (1^9) for the first R-order, the total zonal velocity during the ther- 
mal bar period is 

r, l6jt^A^ ^. 1 ,^ 2 , ^ 1..^ z sinh az . , 
"" -^^1 — 81 — 21? ^ 2az - ~ z smh az - — ^) sm 2itx - sin 1+jrx 

. /2 cosh az - Ix . ^ ^ . cosh 2azv . ^ ^ i/o 
+ (- z smh az + — ^ ) sin birx + (1 - -^ ) sin Sttx] + £i / 

• [b (X) /2 e-fe^2 sin ^ ^^^2^V ^-/^^ (eos ^ . sin ^) - R vA2b,(x)e-^^2 

° \f2 42 42 42 

• sin — ] -eV-=[ (b (Tii)e'''-^ sin ^ + b (112)6 ^^' sin -^) 

42 4k ° 42 " 42 
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- R — (b3(Tii)e^^^2 sin -^ - t3(Tte)e"^^^^ sin -^] + oCl^eV^) (178) 
nTk -^2 42 

where b (x), bi(x), b (t]), b3(T)) have shown in equations (114), (I58), (l2l), 

(IU7). The total vertical velocity is, 

R52rt'*A^ r/ 2(z cosh az - G) 2(sinh az - H) sinh az - Hg z-J-s 
w = ei{ -Q^^ [( - - :^ - ^-^^ + -^^) 

, ,z cosh az - G sinh az - H^ ^ , /2(z cosh az - G) 

' cos 2roc + ( - -—5 ) 2 cos 4jtx - ( 

Ha Ha"^ Ha 

2(shin az - H) cosh az - G z - 1. , ,sinh 2az - Ho , 

- — ^ 5 + 5 - 3-) 5 cos 6itx + ( 5 ^ -z + 1) 

Ha^ 2H^a 2H^ ^ "^ ^ 2H^a '^ 

• i+ cos 8nx] + b '(x) [e"^^^^ cos -^ - 1 + e"^^^ (cos -^ + sin -^)] 

■/2 n/'2 42 

- Rbi'(x) e-^'^^ (cos -^ + sin i£)}-£V2^ [(b (Tii)e^^^^ sin -^ 

n/"2 "/2 4Y.y ° 'v/'2 

nA2 n/2 ^2 

+ o(R2eV2) (1^9) 



where 



8jt^A 
b'(x) = (cos 2jtx - 2 cos hnx) (180) 



and 



, , , ^ 52Jt^A^ r.2G 2 Hp ^ 1 , , G 1 ^ 

bi(x) = — op- |^(— _ -^ _ ^-1- + — ^) cos 2jtx + 2 (— - 72) cos te 



i^H'='a 2H''' ^Ha a' 



,2G _ _2_ ^ G-1 _1^ 
^Ha ~ a^ 2H^a " 2E 



H 1 

5) 5 cos 6jtx + (rrf- - 1) 4 cos 8jtx . (181) 
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The total nondimensional zonal and vertical velocities are plotted in Figure 
20. There are four cells circulating in the cross -section of the lake. Two 
large circulation cells between the thermal bars show a general broad upwelling 
in the middle portion of the lake and divergent flow from the center of the 
lake toward the thermal bar on each side of the lake. Two other smaller circu- 
lation cells are located between the thermal bars and their nearest shores. 
The four cells circulate clockwise and counterclockwise in such a way that 
sinking motions occur along the isotherms at the thermal bars and upwellings 
are present in the middle portion of the lake as well as at the both eastern 
and western shores. The zonal flows are almost entirely in the boundary re- 
gions. The flows along the slant-side-boundaries are more intense. The mag- 
nitude of zonal flow in the free surface Ekman layer is slightly stronger than 
that of the bottom layer. The vertical velocity apparently shows little change 
with depth. The maximum zonal velocity may reach 1 cm/sec and the maximum 
vertical velocity is in general two orders smaller. 



VIII. COMPARISON OF THE RESULTS FROM THE THEORETICAL 
MODEL STUDY WITH FIELD OBSERVATIONS 



As mentioned in the Introduction, Lake Michigan is one of the well surveyed 
large lakes in the world. Beginning at the end of the last cent\iry, a great 
many scientists have devoted themselves to the understanding of limnological 
phenomena and physical dynamics in the lake. Both current studies and tempera- 
ture measurements in Lake Michigan are quite well documented. However, due to 
the complexity of the hydrodynamics in the lake and the multi-variability of the 
surrounding environment, detailed knowledge of the current circulation in Lake 
Michigan is still far from satisfactory. Since the thermal structure in the 
lake is more stable and the method for temperature measurement is simpler, much 
detailed and reliable field data of Lake Michigan are available. 

Theoretical Results 

In our theoretical study three types of temperature distributions and 
circulation patterns have been found at different periods of time around the 
year. During the winter- cooling period, the lake has already been cooled down 
to a temperature less than k°C and both edges of the lake are still colder than 
the middle which is around the temperature of maximum density. The temperature 
distribution of this period has a profile shown in Figure 15, but the horizontal 
temperature difference is in the negative sense. This means that the isotherms 
are distributed as shown in Figure 15^ but the temperature of each isotherm is 
successively decreasing from 4°C at the middle to the lowest temperature at the 
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edges* The meridional velocity circulation of this period is the type-A cir- 
culation as sho-wn in Figure 17 with northbound flow on the east side and south- 
bound flow on the west side of the lake. The zonal and the vertical velocity 
components as plotted as in Figure l8 which shows that two circulation cells 
converge at the middle of the lake. The winter- cooling period may last for 
quite a long period of time^ from January to late March. 

As the season progresses, the atmospheric temperature on the land has 
become much warmer than that of the lake in late March or early April. The 
spring- heating period has the temperature distribution shown in Figure 15. The 
meridional current during this period should be the type-B circulation with 
southbound flow on the east side and northbound flow on the west side of the 
lake as opposed to that shown in Figure 17 • Since the middle portion of the 
lake is colder, being at temperatures less than 4°C, than the temperatures of 
edges at nearly k°Cf the zonal and vertical circulation is Just opposite in 
direction to that shown in Figure l8. Note that the magnitudes of all current 
components are small* 

The spring thermal bar period may start in late April or May. The tempera- 
ture distribution is shown in Figure l6 with the middle portion of the lake at 
a temperature still below ^°C. The meridional current is more complicated as 
shown in Figure 19. The four-celled zonal and vertical circulation pattern 
unique to the thermal bar period is shown in Figure 20. 

As soon as the thermal bar disappears, the slimmer- heating period exists 
before the lake becomes strongly stratified. The temperature of the summer- 
heating period is shown in Figure 15 and the currents belong to the type-A 
circulation. 
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Not until the long (nearly five months) strong- stratification period is 
over will the lake return to almost homogeneous conditions. The autumn- cooling 
period may start in late November or early December. The temperature distribu- 
tion of the autumn- cooling period is exactly the same as in the winter- cooling 
period except the temperature of the whole lake is still at a temperature above 
k^'C and the maximum temperature at the middle is much higher than k°C^ The 
currents are typical type-B circulation. 

The autumn thermal bar current^ before the starting of the next cycle^ is 
usually very weak and may possibly be hardly detected. 

Closing the South End of the Lake Michigan Model 

In our modeling considerations^ since Lake Michigan is connected to Lake 
Huron, the length of the model is considered much longer than the width and 
the flow is treated as meridionally independent. Actually the south end of the 
lake is closed and the Straits of Mackinac also imposes some restrictions on 
the flow at the north end of the lake. However, we will notice that the domi- 
nant meridional current is anti- symmetrical to the center line of the lake. 
The meridional currents on each side of the lake are always flowing in opposite 
directions. Once the ends are closed or restricted, the theoretical flows on 
each side of the lake can be smoothly connected together at the ends by the 
continuity requirements^ while the meridional currents become cyclonic or anti- 
cyclonic flows in the lake. We also notice in Figure 2 the existence of an 
everywhere positive north-south meridional temperature gradient. From the re- 



139 

suits obtained above, we deduce that the cyclonic flow in the lake is strengthened 
by this meridional temperature gradient, and the anticyclonic flow is weakened. 

Comparison of Temperatures 

Church {l9h2B.,'b) made a comprehensive study of the annual temperature cycle 
of Lake Michigan in many cross- sections of the lake during 1941-^2. His tempera- 
ture data, to my knowledge, are still the most complete for Lake Michigan. Some 
of his published figures are reproduced here for comparison with our results ob- 
tained from the theoretical study during the corresponding periods of the year. 
Figure 21 is the reproduction of Church's January temperature data across the 
Milwaukee-Muskegon section near the central part of Lake Michigan. As we 
noticed from Figure 21, the maximum surface temperature in the cross-section is 
at i|-°C in the middle of the lake, while curved and concave isotherms decrease 
toward the shores. This is exactly the temperature distribution obtained in 
our theoretical study during the winter- cooling period discussed above, except 
that the former has sharper temperature gradient near the shores. Qualitatively, 
Figure 21 agrees very well with Figure 15. Figures 22 and 25 are the reproduc- 
tions of Churches February and March temperature. They are all of the same 
distribution pattern as given by the theoretical study. 

In late March or early April, the edges of the lake have already been 
warmed up to a temperature nearly ^"^C or slightly above as shown in Figure 2k^ 
which is Church's observed temperature distribution of middle April. Figure 2k 
shows the characteristic temperature distribution of the spring-heating period 
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which has two edges at nearly 4°C and the middle of the lake at slightly lower 
temperature. Qualitatively the theoretical prediction agrees well with the 
field data. 

As the atmospheric temperature gets higher and higher, the thermal bar 
period immediately follows. Figure 2k already shows some of the "beginning 
phenomenon of a thermal "bar at each edge. Detailed temperature distribution 
across a thermal har is not available from Churches temperature figures. 

Rodgers (1965^ I966) has investigated the thermal "bar in the Great Lakes 

and has measured the temperature across the thermal "bar in great detail. His 

temperature section across the thermal "bar is reproduced here as shown in 

Figure 25. Notice from Figure 25 the temperature near the east shore is above 
WEST EAST 
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Figure 25 » Temperature distribution across a section of 
the thermal bar (after Rodgers, I965). 
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12°C and the temperature of the open lake indicated by the "west" end in the 
figure is less than the temperature of maximuin density. The temperature between 
two broken isotherms shown in the figure is the band of fluid of nearly ^°C 
called the thermal bar- There is a strong horizontal temperature gradient 
toward the shore and a weak one toward the center of the lake. Isotherms be- 
tween the thermal bar and the shore are concave up and shoreward. Noble and 
Anderson (I968) found the same temperature pattern across a thermal bar off 
Grand Haven in Lake Michigan as that of Rodger s in Lake Huron. The temperature 
distribution across a thermal bar shows exactly the features of our theoretical 
temperature solution duiring the thermal bar period obtained from the model study 
as shown in Figure I6. 

The thermal bar period usually lasts four to six weeks. After the disap- 
pearance of the thermal bar the summer- heating period starts. The surface 
temperature of the lake is above 4*^0 everywhere at this stage with the minimum 
temperature near 4°C at the middle of the lake. Figure 26 is the reproduction 
of Churches May temperature data which shows the general feature of the summer- 
heating period and qualitatively agrees well with the theoretical temperature 
distribution of this period as shown in Figure 15 . Figure 27 is the June 
temperature distribution of a cross section farther north of Lake Michigan^ 
which still shows the basic temperature structure of the summer-heating period 
indicated in Figure I5. 

The thermocline in Lake Michigan usually starts to form not long after the 
disappearance of the thermal bar and gets stronger in June as indicated in 
Figure 27. At the beginning^ the thermocline may be very weak. The lake may 
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form a shallow thermoeline during the hot day which disappears at night* As 
the temperature gets higher and higher^ the thermoeline gets stronger and 
stronger after June. Not until the thermoeline is very deep and weak in late 
November as indicated in Figure 28 ean the lake be considered in a homogeneous 
state again. At that time the atmospheric temperature has already become colder 
than the temperature of the lake. 

After the disappearance of the thermoeline^ the autumn- cooling period be- 
gins. The temperature at the middle portion of the lake is higher than that of 
the edges as shown in Figure 29. Qualitatively, Figure 29 also shows the pre- 
dicted temperature feature indicated in Figure 15 with temperature decreasing 
toward shores. 

There are no data available to compare the temperature distribution during 
autumn thermal bar period. It will happen theoretically after the autumn- 
cooling period when the edge temperature of the lake is already below ^''C. As 
we have mentioned, due to the smallness of the horizontal temperature differ- 
ences and to the instability of the air mass above the lake the autumn thermal 
bar may be hardly observed. 

In general, the predicted temperature distributions of our model study 
agree qualitatively very well with the observed field data. But the isotherms 
of the theoretical temperature solution are more widely spread and evenly dis- 
tributed than the real data. The location of the thermal bar is also farther 
from shore than observed. This is because of the smoothed surface temperature 
profile we used as a boundary condition in solving the energy equation. 
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Comparigon of Velocity Field 



MERIDIONAL VELOCITY 

The configuration of surface current during the navigation season found by- 
Harrington (1895) by means of drift bottles (drawn by Millar^ I952) was repro- 
duced in Figure 1. His general mean current pattern in Lake Michigan shows a 
cyclonic circulation. There are also two cyclonic cells together with the 
general lake circulation^ one well developed cell in the southern basin and 
another weaker cell in the northern basin. The northward current on the east 
side of the lake is very well defined. The southward current on the west side 
of the lake is weaker. The current in Lake Michigan varies in speed from 10 to 
20 cm/ sec. 

Harrington's current pattern has no specific time correspondence with our 
defined time periods. However^ the summer-heating period^ the spring thermal 
bar period^ the strong stratification period^ and the autumn- cooling period are 
within the navigational season. Of course^ the strong stratification period is 
out of our prediction. Of the other three periods^ the two longer periods show 
cyclonic circulation in the lake. Duri^ig the summer- heating period^ the meridi- 
onal velocity shpwn in Figure 17 gives a surface current pattern like that of 
Harrington. During the thermal bar period^ the meridional velocity is shown in 
Figure I9 which shows that the currents close to both shores have cyclonic cir- 
culations while the middle portion of the lake has anticyclonic circulation. 
The circulation of the autumn- cooling period is reversed. In our model study^ 
the current magnitudes may reach a maxim.ura of 20 cm/sec during the summer- 



152 

heating season^ 12 cm/sec during the spring thermal "bar period which coincide 
with Harrington's measurements. However, due to the lack of precise time 
correspondences, the comparison between our theoretical solutions and Harrington's 
field results is very loose. The meridional circulations of the spring- thermal 
bar and the summer- heating periods agree qualitatively with Harrington's, but 
the autumn- cooling period does not. Furthermore, our simple model study does 
not show detailed eddy patterns in the lake. 

Ayers et al. (1958) made four cruises pf synoptic survey in Lake Michigan 
in 1955 and studied the current by using drift bottles, then computed the cur- 
rent circulation by the modified dynamic height method (Ayers, 1957) • Their 
August current circulation, while the lake is strongly stratified at that time, 
is not predicted in our homogeneous model study. Their June surface currents 
have been reproduced here as shown in Figure 50* Though a relatively weak 
thermocline has already been formed during the time of their cruises, our 
summer-heating period circulation pattern should still be observed if there is 
no drastic change in the lake dynamic system. Ayers' June currents show, in 
general, a cyclonic circulation with two weaker cells in the northern and 
southern basin similar to that of Harrington's; but in addition, they found 
another large (clockwise gyre near Benton Harbor, Michigan City area at the 
southern Lake Michigan basin. The magnitudes of currents they found are generally 
less than 10 cm/sec. Our meridional current of the surmner- heating period is 
shown in Figure 17. Qualitatively, the general circulation pattern and speed 
magnitude agree well with Ayers' surface currents. Note that the detailed 
eddy circulations have been neglected in our model study. 
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Figure 50. June sTirface currents in Lake Michigan (from Ayers et al., I958). 



15^ 

In 1962 and 196^, the FWPCA deployed numerous temperature and current re- 
cording instruments at more than thirty stations all over the lake. At each 
station there was an instrument pair consisting of temperature recorder and 
current meter at each fixed depth of 10^ 15^ 22, 50 m and at each succeeding 
30-m level. Unfortunately, these finely spaced, long period, lakewide recorded 
data have not been analyzed and published. Their preliminary report (U. S. 
Dept. of Interior, I967) indicates that, generally speaking, surface currents 
in Lake Michigan are weak and southbound near the west shore, more intense and 
northbound near the east shore. No permanent type of circulation pattern has 
been found in Lake Michigan, but four basic patterns have been observed, two in 
winter and two in summer. Neglecting the summer patterns, their roughly defined 
winter surface circulation patterns are reproduced as shown in Figures 5I ^^^ 
52. Figure 5I is the winter pattern which normally occurs from January to April 
though not necessarily in a continuous manner. This pattern accounts for 20 to 
25 percent of the water movements during the year. This current pattern shows 
a general cyclonic circulation with an elongated gyre on the northern basin 
and another large gyre in the southern basin. The corresponding time in our 
categorization is the winter- cooling period which has shown in general the same 
circulation pattern indicated in Figure I7 as that in Figure 5I. The other 
winter pattern they reported is shown in Figure 52 which is usually found after 
November through March intermittently. This time period is relevant for compari- 
son with our autumn- cooling period which occurs after the thermocline is very 
deep and weak, or after it disappears. They both show an anticyclonic circu- 
lation. The eddy gyres and nearshore flows of the winter patterns do not agree 
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Figure 31. First Tmsl^ -winter stirfaee current pattern in 
lake Mlelilgan (tsmm U* S* lept^ #f Interltrji 196t)» 
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Figure 52* Second "basic winter surface eurrent pattern In 
Lake MicM^n (from II, S* Deptr of Interior, I96T). 
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well with our theoretical prediction. This may be due to our neglect of wind 
stresses which is important to fluctuating currents^ especially in the very 
shallow area as well as the uppermost surface layer of the lake. 

On 10 May and 2J June 1968^ we ran two sets of drogues for the current 
study in the west side of Lake Michigan outside of Milwaukee Harbor as indicated 
in Figure 55. The wind conditions^ positions of each fixed pointy running time ^ 
and current speed are shown in Table I5. Currents were flowing southward in 




43^*00 



MilwouliM 



CALM 



E. 



(a) 



/ -. 



87**A0' 

F. 



.43 



• i 



87*50' 



43**00» 



MiUt 



T 





(b) 



Figure 53. Drogue studies for swctace current measurements near western 
shore. (a) 10 May I968. (b) 27 June I968. 
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TABLE 15 
Drogue Studies of Surface Current on the West Side of Lake Michigan 



Date 


10 May 1968 


27 June 1968 


Wind 


N, k kt 


NxW, 6 kt 


Lake 
Condition 


Ca.1 m 




1 ft Wave 




Station 


Time 


Speed, 
cm/ sec 


Mean, 
cm/ sec 


Station - 


Time 


Speed, 
cm/ sec 


Mean, 




Start 


Finish 


Start 


Finish 


cm/ sec 




A 


0800 


1506 


5.12 




H 


1245 


1848 


4.91 






B 


0815 


1520 


iv.75 




I 


2151 


1820 


9.86 




Drogue 
Study- 


C 


0821 


1552 


8.02 




J 


1301 


1813 


10.74 




D 


0837 


l^kQ 


i^.91 


5.59 


K 


1308 


1805 


7.14 


7.6 




E 


0848 


l66k 


7.58 




L 


1314 


1856 


4.91 






F 


0859 


1620 


5.56 




M 


1321 


1840. 


8.47 






G 


0912 


1655 


4.01 




N 


1327 


1831 


7.14 





agreement with our meridional currents on the west side of the lake during the 
summer- heating period. On the east side of the lake^ a tripod- supported pendu- 
lum current meter was installed 250 m offshore at a depth of 5 m from 11 May to 
11 October I967 under Dr. Ayers' direction. The late spring and early summer 
current directions are summarized in Table ik. We notice that most of the cur- 
rents flow northward near the east shore of the lake during the summer- heating 
period, which agrees well with our meridional velocity solution prediction. 

TABLE Ik 



Percentage of Hours Current Flowed Northward^ Westward, and Southward on the 
East Side of Lake Michigan Near Benton Harbor, Michigan 



Month 




Current , ^ 




Northward 


Westward 


Southward 


]yfa,y 


96.7 


5.0 


0.5 


June 


95.2 


k.& 





July 


97.1 


2.9 






1^9 

ZONAL AND VERTICAL VELOCITIES 

Due to the smaliness in magnitudes of the zonal and vertical velocities 
there are not much field data availabij-e for comparisons with our theoretical 
current circulations in th^ cross-section of Lake Michigan. However^ the thermal 
bar phenomenon in a large dimictic lake is quite well observed and documented 
by Tikhomirov (1965)^ Rodgei^i? (1965^ 1966)^ and others. The thermal bar is 
established in the spring-warming-^p seaspn when the shallow inshore edges of 
the lake have been warmed to a]?Qve k^C, while the deep open- lake-water is still 
less than k^'C. At the junction of the two waters, a boundary band of water (the 
thermal bar) of nearly k^'C is created by mixing of the warmer inshore water and 
the colder offshore water. The 4°C water fQ?:*med at the thermal bar is denser 
and it sinks. The sinking water is replaced by a shoreward surface flow from 
the deep open- lake toward the thermal bar^ and by a lakeward surface flow from 
the inshore region toward the thermal bar. The convergence zone is readily ob- 
served by the windrows of dead fish, trash, and garbage accumulated along the 
thermal bar. 

Rodger s (I965) has deduced a thermsil bar circulation as shown in Figure 
5^. There are four celjs circulating contratotationally in the cross section 
with sinking motion alpng the ^''C isotherms at the thermal bars and upwellings 
at the middle as well as on both edges of the lake. This kind of current pat- 
tern is exactly our thermal bar circulation as shown in Figure 20. Note that 
we use a smoother temperature profile in our model investigation which results 
in moving the thermal bar farther ovit to the open lake. 
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Figure jU. Thermal bar circulation (after Rodger s, I965). 



The actual condition of the thermal bar mixing phenomenon is considerably- 
more complex as indicated in the aerial photograph (Figure 55) taken along the 
thermal bar near Grand Haven in Lake Michigan during the spring of I968. Zonal 
and vertical current meastirements of the thermal bar circulation are very dif- 
ficult, not only because of their small magnitudes but also because of their 
being obscured by other factors. 
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IX, CONCLUSIONS 

We have theoretically Investigated the thermally induced current in Lake 
Michigan. Using the procedure of douhle perturbation expansions outlined in 
Chapter III^ approximate solutions have "been obtained for the vhole flow field, 
except for four unimportant corners • Approximation errors of the solutions 
are estimated to "be the order of (Rosshy number) . 

Two types of nondimensional temperature distributions have been found. 
One type prevails during all the general thermal current periods, namely winter- 
cooling, spring-heating^ summer-heating, and autumn-cooling periods, with the 
temperature distributions as shown in Figure 15* The other type prevails dur*- 
ing thermal bar periods, both spring and autumn^as shown in Figure l6. Three 
types of circulation pattern have been found in the model study. There are 
the type-A and the type-B circulation patterns^ the general thermal currents 
(Figures 17 and l8), and the type-C circulation pattern (Figures I9 and 20), 
The type-B currents are exactly the opposite of the type-A currents. The type- 
C circulation is unique for the thermal bar phenomena. 

In trying to verify the theoretical thermal circulation currents, some 
field measurements for currents and eddy viscosities were conducted in I967 
and 1968 in Lake Michigan. Comparisons have been made between predictions ob- 
tained from our mathematical study and previously published field data such 
as those of Harrington (l895)> of Church (igii-ga^b)^ of Ayers et al« (1958)^ 
and of Rodgers (I965)* Qualitatively, the theoretical temperature distribu- 
tion and the thermal current structure agree very well with the real field 
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olDservation^ in lake Michigan during the corresponding seasons. This leads 
one to postulate th^^t the thermal hody force does play an important role in 
generating th? mean lake circulation patterns. 

Our theoretical solutions do not show large eddies in the southern and 
northern "basinig which have heen observed. This is perhaps heoause the model 
assumes a flat "bottom* 

The agreements "between the theoretical solutions and the field ohserva*- 
tions are rather poor in the shallow regions^ This may well he the result of 
the faster response of the shallow water current to wind stresses and of the 
irregularities of the coastlines., Due to the neglect of wind stresses^ the 
model current circulation is not expected to coincide with instant current 
measurem,ents> especially in the shallow water or at the upmost surface layer^ 
hut it is expected to he valid for the prediction of the mean current cir^ 
culation pattern in the lake. 

As to the thermal currents in lake Michigan during the strongly strati- 
fied period in the summer and early autu^n^ further theoretical work needs 
to he done. 



REFERENCES 



Ayers^ J« G. 1957* A dynamic height method for the determination of currents 
in deep lakes. Limnology and Oceanography _, 2:150-l6r. 

. 1962. Great Lakes Water^ their circulation and physical and 



chemical characteristics. Amer. Assoc. Adv. Sc*> Puhl, No Tl: 71-89- 
, D. C. Chandler^ G. H. Lauff^ and C, F. Powers. 1958* Currents 



and water masses of Lake Michigan* Univ, of Michigan^ Great Lakes Re- 
search Division^ Puhl, No. 5^ I69 p. 

Barcilon^ Victor. 196^^ Role of the Ekman Layers in the stahility of the 

symmetric regime obtained in a rotating annulus. J. of the Atmospheric 
Sciences^ 21:291-299. 

Batchelor^ G. K. 1953 • The theory of homogeneous turbulence^ Cambridge 
University Press^ I95 p* 

and A. Ao Townsend. I956. Turbulent diffusion. Surveys in 



Mechanics c Cambridge University Press^ pp. 552-599. 

Bellaire^ F. R. 19^5 . The modification of warm air moving over cold water. 
Proc^ 8th Conf . on Great Lakes Research^ Univ,. of Michigan^ Great Lakes 
Research Division Publ. No. 15^ pp. 2^9-256. 

Boussinesq.^ J« I877. Theorie de l*ecoulement tourbillant. Mem. Proc. Acad. 
Sci, XXII> k6. Paris^ 

Bowden^ K. F. 196^* Turbulence* Oceanographic and Marine Biology Annual 
Review^ 2:11. 

j> K« F« 1965. Horizontal mixing in the sea due to a shearing current. 



Jo Fluid Mech.^ 21:85-95. 

Bowles^ P*^ R* H. Burns^ F* Hudswell^ and R.R.P* Whipple* I958. Sea dis- 
posal of low activity effluent* Proe. 2nd International Conf. Peaceful 
Use of Atomic Energy* United Nations^ N* Y.^ 18:576-589. 

Brindley^ J. 1960^ Stability of flow in a rotating viscous imcompressible 

fluid subjected to differential heatingo Royal Soc* London^ Phil. Trans. ^ 
A255:l-25. 

Carrier^ G. T. 1955« Boundary layer problems in applied mechanics, Adv» in 
Applied Mech*^ 5:1-19. 

I6i^ 



X65 



Charnock, H. X959. Tida.! friction from currents near the sea "bed. Geophys. 
J., 2:215-221- 

Churchy P, E. 19^2a, ThQ annual temperature cycle of Lake Michigan. Univ* 
of Chicago^ Institute of Meteorology^ Misc. Report No. k^ kQ p. 

19J^2"b, The annual temperature cycle of Lake Michigan. Univ* of 



Chicago^ Institute of Meteorology^ Misc. Report No- l8^ 100 p, 

Csanady^ G^ T, I963, Turbulent diffusion in Lake Huron* J* Fluid Mech*^ 
17:560-585* 

196^, Turbulence and diffusion in the Great Lakes, Univ, of 



Michigan^ Great Lakes Research Division Pahl. No, 11:526-559, 
. 1966, Accelerated diffusion in the skeved shear flow of lake 



currents. J. Geophys. Res.^ 71:^ll-^20» 
. 1967* Large-scale motion in the Great Lakes^ J. Geophys^ Res,^ 



72:^151-^162. 



, 1968a. Motion^ in a model great lake due to a suddenly imposed 



wind. J. Geophys. Res.^ 75: 6^1-55 -6Uii-7* 
, 1968h* Wind -driven summer circulation in the Great Lakes. J, 



Geophys. Res*^ 75:2579-2589* 

Davies^ T. V. 1956* The forced flpw due to heating of a rotating fluid. 
Royal Soc, London^ Phil* Tran3.^ A249:27-6i^.. 

^ 1959* On the forced motion due to heating of a deep rotating 



liquid in an annulus. J. Fluid Mech.^ 5:595-621* 

Davis^ R. A, and D^F^R, McGeary. 1965, Stability in near-shore bottom topog- 
raphy and sediment distribution, southeastern Lake Michigan, Proc. 8th 
Conf * on Great Lakes Research^ Univ* of Michigan, Great Lakes Research 
Division Publ<. No- 15:222-251. 

Deacon, E* L. 1962* Aerodynamic roughness of the sea, J. of Geophys* Res,, 
67:5157-5172, 

Deason, H» J* 1952. A study of surface currents in Lake Michigan* The 
Fisherman, 1:5--^, 12. 

Defant^ A* I961, Physical oceanography* MacMillan Co*, N, Y#, Vol* I, 729 



166 



Eady^ E* T. 19^9* Long waves and cyclone waves. Tellus^ 1:33-52. 

Ekman, V. W, I905. On the influence of the earth ^s rotation on ocean currents* 
Arkiv. Mat^ Astr. Fysik,;, 2:1-52* 

Elder^ F. C* and E. S. Soo, I967* An investigation of atmospheric turhulent 
transfer processes over water. Univ, of Michigan^ Great Lakes Division 
Special Report No. 29^ k6 p. 

Farlow^ J. S. I965. A field technique used for horizontal diffusion studies 
in Lake Michigan and Lake Erie* Proc, 8th Conf. on Great Lakes Research^ 
Univ. of Michigan^ Great Lakes Research Division Pahl. No. 13:299*'503« 

Fowlis^ W« ¥^ 1963. An experimental study of the transitions between the 

flow regime^ of thermal convection in rotating annulus of liquid* Ph»D« 
Thesis^ Univ* of Durham^ England^ I88 p. 

and R. Hide. I965. Thermal convection in a rotating annulus of 



liquid; effect of viscosity on the transition between axisymmetric and 
non-symmetric flow regimes, J. Atmosph. Sc.^ 22:5^1-558. 

Foxworthy^ J. E. and G. M. Barsom. I967. Comments on paper by G. T. Csanady 
"Accelerated diffusion in the skewed shear flow of lake currents. " J. 
Geophys. Res., 72:2683-2691. 

Fultz, D. 1953 » A survey of certain thermally and mechanically driven fluid 
systems of meteorological interest. Fluid models in geophysics, Proc. 
1st Symposium on the Use of Models in Geophysical Fluid Dynamics. Balti- 
more, Maryland, pp. 2T-63, 

Greenspan, H. P. and L. N. Howard. I963. On the time -dependent motion of a 
rotating fluid. J^ Fluid Mech., 17:385-^01^-. 

Hadley, G. 1735. Concerning the cause of the general trade winds. Royal 
Soc. London, Phil. Tans., 59:58-63. 

Harrington, M. W. I895. Surface currents of the Great Lakes. U. S. Dept. of 
Agriculture, Weather Bureau, Bull* B. 

Hide, R. 1953. Some experiments on thermal convection in a rotating liquid* 
Royal Meteorological Soc, Quart. Jour., 79:l6l. 

1958. An experimental study of thermal convection in a rotating 



liquid. Royal Soc. London, Phil. Trans., A250: ^^1-^^-78. 
. 196^4-. The viscous boundary layer at the free surface of a rotating 



baroclinic fluid. Tellus, 16:523-529- 



167 



Hill, M. N. Ed. 1962. The sea. Vol. I. Physical oceanography. Interscience, 
Wiley, N. Y. , 86i|- p. 

Hinze, J. 0. 1959- Turhulence. McGraw-Hill, N. Y., 586 p. 

Hough, J. L. 1958. Geology of the Great Lakes. University of Illinois 
Press, Urhana, 111., 313 P« 

Hunt, I. A. 1959. Winds, wind set-ups, and seiches on Lake Erie. U. S. Lake 
Survey, Detroit, Michigan, 59 P* 

Hunter, C. 196?. The axi symmetric flow in a rotating annulus due to a hori- 
zontally applied temperature gradient. J. Fluid Mech., 27:753-778. 

Hutchinson, G. E. 1957. A treatise on limnology. Vol. I. Geography, 
Physics, and Chemistry. Wiley, N. Y., 1015 p. 

Kolmogoroff, A. N, 19^1. The local structure of turhulence in an incompres- 
sihle viscous fluid for very large Reynolds numher. Doklady Akad. Nauk 
U.S.S.R. Geophys., Ser. 3: 1+05-^13. 

Kuo, H. L. 1957. Further studies of thermal motions in a rotating fluid. 
J. Meteorology, 1^:553-558. 

Lamb, H. 1932. Hydrodynamics. Dover, N^ Y. , 673 p. 

Lesser, R. M. 1951. Some observations of the velocity profile near the sea 
floor. Amer. Geophys. Union, Trans. 32:207-211. 

McFadden, J. D. and R. A. Ragotzkie. 1963. Aerial mapping of surface temper- 
ature pattern of Lake Michigan. Proc . 6th Conf. on Great Lakes Research, 
Univ. of Michigan, Great Lakes Research Division Puhl. No. 10:55-58. 

Millar, F. G. 1952- Surface temperatures of the Great Lakes. Fish Res. Bd. 
Canada Jour., 9:329-376. 

Mortimer, C. H. I963. Frontiers in physical limnology with particular ref- 
erence to long waves in rotating "basins. Proc. 6th Conf. on Great Lakes 
Research, Univ. of Michigan, Great Lakes Research Division Puhl., 10: 
9-^2. 

Munk, W. H. 1950. On the wind-driven ocean circulation. J. Meteorology, 
7:79. 

Noble, V. E. 1966. Vertical current structure in the Great Lakes. Univer- 
sity of Michigan, Great Lakes Research Division Special Report No. 27, 
i^2 p. 



168 



. 1967. Evidences of geostrophically defined circulation in Lake 

Michigan, Proc. 10th Conf. on Great Lakes Research, pp. 289-298. Inter- 
national Association for Great Lakes Research. 

Noble, V. E. and R. F. Anderson. I968. Temperature and current in the Grand 
Haven, Michigan vicinity during thermal "bar conditions. Proc. 11th Conf. 
on Great Lakes Research, pp. U70-U79. International Association for 
Great Lakes Research. 

Noble, V. E. and J. C. Ayers. I96I. A portable photocell fluorometer for 
dilution measurements in natural waters. Limnology and Oceanography, 
6:I^57-i+6l. 

Okubo, A. 1962, A review of theoretical models for turbulent diffusion in 
the sea. J. Oceanographical Society of Japan, the 20th annual volume, 
286 p. 

and S. Farlow. 1967* Analysis of some Great Lakes drogue stud- 



ies. Proc, 10th Conf on Great Lakes Research, pp. 299-508. Interna- 
tional Association for Great Lakes Research. 

Pierson, W.J., Jr. 196^. The interpretation of wave spectrum in terms of 
the wind profile instead of the wind measured at a constant height. 
J. Geophys. Res., 69. 

Richardson L. E. I926. Atmospheric diffusion shown on a distance neighbor 
graph. Proc. Royal Soc. London, Phil. Trans., A110:709-727. 

and H. Stommel. 19^8. Note on eddy dlffusivity in the sea. J. 



Meteorology, 5:258-2^^0. 

Richardson, W. S. I962. Instruction manual for recording current meter. 
Woods Hole Oceanographic Instituttion, Ref. No 62-6, 10 p. 

Robinson, A. R. 1959* The symmetric state of a rotating fluid differentially 
heated in the horizontal. J. Fluid Mech., 6:599-620. 

Rodgers, G. K. 1965* The thermal bar in the Laurentian Great Lakes. Proc 
8th Conf. on Great Lakes Research, Univ. of Michigan, Great Lakes Re- 
search Division Publ. No. 15: 558-565 • 

. 1966. The thermal bar in Lake Ontario, spring 1965^ and winter 



1965-66. Proc. 9th Conf, on Great Lakes Research, Univ. of Michigan, 
Great Lakes Research Division Publ. No. . 15:569-57^. 

Schllchting, H. I968. Boundary- layer theory. McGraw-Hill, N. Y. , 757 P. 



X69 



Stoermer^ E. F* I968, Near shore phytoplankton populations in the Grand Haven^ 
Michigan vicinity during thermal har conditions, Proc* 11th Conf . on 
Great lakes Research^ pp* I37-I5O. International Association for Great 
Lakes Research* 

Storamel^ H, 1955. Lateral eddy viscosity in the Gulf Stream system. Deep- 
Sea Res.> 5:88-90. 

Strong, A, E. and F. R. Bellaire. 1965. The effect of air stability on wind 
and waves* Proc. 8th Conf. on Great Lakes Research^ Univ. of Michigan^ 
Great Lak^s Research Division Puhl. No, 13:285-289. 

Sverdrup^ H* U«, M* W. Johnson, and R. H. Fleming. 19^4-2. The oceans. Their 
physicsj chemistry, and general biology. Prentice-Hall^ N, Y. , IO87 p, 

Tikhomirov, A* I. I965. The thermal bar of Lake Ladoga^ Bull* (izvestiya) 
All -Union Geograph, Soc. 95> American Geophys. Union Transit Soviet 
Hydrology^ selected paper No. 2:182-191* 

U. S. Department of Commerce, Weather Bureau. 1959* Climatology and weather 
services of the St. Lawrence Seaway and Qreat Lakes* Prepared by Marine 
Area Section, Office of Climatology, U, S. Dept. of Commerce, Weather 
Bureau Teclinlcal Paper No* 55 > 75 P* 

U. S. Department of the Interior, Federal Water Pollution Control Adminis- 
tration, 1967. Lake current s# A technical report containing back- 
ground data for a water pollution control program. Water Quality In- 
vestigations, Lake Michigan Basin, U. S. Dept. of the Interior, Federal 
Water Pollution Control Administration, Great Lakes Region, Chicago, 
111., 564 pp. 

Van Dyke, Milton. I96U. Perturbation methods in fluid mechanics. Academic 
Press, N. Y.^ 229 p. 

Verber, J* L. 196^^-. Initial current studies in Lake Michigan. Limnology 
and Oceanography, g:k26'^k'}0^ 

* 1965. Current profiles to depth in Lake Michigan. Proc* 8th 



Conf. on Great Lakes Research, Univ. of Michigan, Great Lakes Research 
Division Publ. No. 15:564-571. 

Welch, P. S. 1952. Limnology. 2nd ed. McGraw-Hill^ N» Y., 558 p. 



